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The  potential  flow  about  a  prolate  spheroid  in  axial  horizontal 

* 

motion  beneath  a  free  surface  is  treated  analytically.  While  the  free- 
surface  boundary  condition  is  linearized,  the  boundary  condition  on  the 
surface  of  the  body  is  satisfied  exactly.  Thus,  an  "exact"  solution, 
within  the  theory  of  infinitesimal  waves,  is  obtained.  The  solution 
is  sought  in  the  form  of  a  distribution  of  sources  on  the  surface  of 
the  spheroid,  of  unknown  density;  the  analysis  yields  an  infinite  set 

I 

of  equations  for  determining  the  coefficients  of  the  expansion  of  the 
density  function  in  spherical  harmonics  (and  therefore  for  determining 
the  coefficients  of  the  expansion  of  the  potential  in  spheroidal  har¬ 
monics).  An  expression  is  derived  for  the  wave  resistance  of  the 
spheroid  in  terms  of  these  coefficients  through  application  of  the 
Lagally  theorem.  The  expression  for  the  wave  resistance  given  by 
Havelock  in  1931  is  obtained  as  the  first  approximation  in  the  present 
analys is . 

Numerical  evaluations  were  performed  using  an  IBM  350/65  computer, 
for  a  Froude  number  (defined  with  respect  to  the  distance  between  foci) 
of  0.4,  a  focal  distance  equal  to  twice  the  depth  of  submergence,  and 
several  values  of  the  eccentricity.  The  numerical  solution  of  the  in¬ 
finite  set  of  equations  was  obtained  keeping  an  increasing  number  pf 
equations  (and,  therefore,  calculating  an  increasing  number  of  coef¬ 
ficients  of  the  series  expansions),  up  to  a  maximum  of  twenty.  The 
wave  resistance  and  t^e  density  of  the  source  distribution  were  eval¬ 
uated  at  each  stage,  the  latter  along  meridian  planes  of  the  spheroid. 

For  a  prolate  spheroid  with  a  slenderness  ratio  slightly  larger  tha  . 

j 

five  the  wave  resistance  is  larger  than  Havelock's  by  about  347,.  For 
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slenderness  ratios  of  3.64  and  2.40  the  corrections  are  as  much  as  907. 
and  368%,  respectively,  of  Havelock's  approximation  (the  spheroid  corres 
ponding  to  the  latter  slenderness  ratio  is  very  close  to  piercing  the 
free  surface). 

An  infinite  set  cf  equations,  essentially  equivalent  to  that  ob¬ 
tained  in  this  work  for  determining  a  source  distribution  on  the  sur¬ 
face  of  the  spheroid  which  satisfies  exactly  the  boundary  condition  on 
its  surface,  was  obtained  by  Bessho  using  an  entiraly  independent  deri¬ 
vation.  The  coefficients  of  Bessho's  system  of  equations,  however, 
appear  to  be  incorrect,  possibly  because  of  typographical  errors,  and 
his  numerical  evaluations  are  rather  inaccurate.  The  value  of  the 
wave  resistance  obtained  by  Bessho,  for  a  Froude  number  of  0.395,  a 
focal  distance  equal  to  twice  the  depth  of  submergence,  and  a  slender¬ 
ness  ratio  of  4.17,  exceeds  Havelock's  approximation  by  1467.;  according 
to  the  numerical  evaluations  reported  here,  the  correction  should  in¬ 
stead  be  in  the  neighborhood  of  60  to  65%  of  Havelock's  value. 
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INTRODUCTION 


Even  the  simplest  problems  involving  free  surface  waves  are 
difficult  to  treat  analytically  when  formulated  exactly.  In  dealing 
with  these  problems,  therefore,  one  must  introduce  simplifying  assump¬ 
tions  into  the  equations  of  motion  and  the  boundary  conditions,  and 
replace  the  original  problem  by  a  new  one  which  is  more  amenable  to 
mathematical  treatment  and  which  should  approximate  the  original  one  in 
some  prescribed  way.  If  viscosity  is  neglected  and  irrotational  flow 
is  assumed,  the  problem  reduces  to  finding  solutions  of  the  Laplace 
equation.  The  free-surface  boundary  condition,  however,  is  nonlinear, 
even  if  surface  tension  is  neglected,  and  moreover  it  must  be  applied 
on  a  surface  of  unknown  location.  As  a  result,  further  simplifica¬ 
tions  are  still  needed  in  most  cases.  For  flows  about  submerged  ob¬ 
stacles,  the  assumption  of  irrotational  flow  may  lead  to  a  new  problem 
which  does  not  approximate  the  original  one;  such  is  the  case,  for  exam¬ 
ple,  with  flow  about  circular  cylinders  and  spheres,  these  two  shapes 

being  very  attractive  to  the  mathematical  researcher  because  of  the 
* 

simplifications  they  bring  about  in  the  mathematical  expressions. 

In  order  to  deal  with  the  nonlinear  free-surface  boundary  condition, 
special  approximation  techniques  have  been  developed.  For  flows  about  ■ 
submerged  bodies,  these  techniques  transform  the  problem  into  a  series 
of  linear  problems  with  nonhomoganeous  boundary  conditions  (with  the 
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exception  of  the  first  order  or  infinitesimal  wave  approximation,  for 
which  the  corresponding  boundary  condition  is  homogeneous),  to  be  ap 
plied  on  the  plane  of  the  undisturbed  free  surface.  Even  within  the 
theory  of  infinitesimal  waves,  however,  an  exact  solution  is  difficult 
to  obtain.  An  additional  approximation  is  involved  in  many  researches, 
since  the  boundary  condition  on  the  surface  of  the  body  is  not  satis¬ 
fied  exactly.  If  one  assumes  that  the  shape  of  the  body  is  such  that 
perfect  fluid  theory  can  yield  results  of  practical  significance,  the 
determination  of  the  errors  due  to  the  failure  to  satisfy  the  nonlinear 
frce-surface  boundary  condition  on  the  one  hand,  and  the  boundary  con¬ 
dition  on  the  surface  of  the  body  on  the  other,  is  of  great  practical 
interest. 

In  the  present  work,  the  wave  motion  produced  when  a  stream  of 
constant  velocity  is  incident  upon  a  submerged  prolate  spheroid  with 
its  axis  parallel  to  the  free  surface  and  in  the  direction  of  the 
stream  will  be  treated.  This  shape  is  of  interest  in  connection  with 
the  wave  resistance  of  submarines;  moreover,  perfect  fluid  theory  can  be 
expected  to  yield  results  of  practical  significance  when  applied  to 
determine  the  flow  about  slender  prolate  spheroids.  The  free  surface 
boundary  condition  will  be  linearized;  the  boundary  condition  on  the 
surface  of  the  body  will  be  satisfied  exactly.  Thus,  an  answer  will 
be  provided,  for  this  particular  shape,  to  the  latter  of  the  two  pro¬ 
blems  posed  in  the  preceding  paragraph. 

A  short  summary  of  past  researches  dealing  with  the  theory  of 
infinitesimal  waves  in  flows  about  submerged  obstacles  is  included  be¬ 
low.  Only  streams  of  infinite  depth  are  considered. 
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The  first  detailed  evaluation  of  the  wave  motion  caused  by  a  fully 
submerged  obstacle  is  contained  in  a  paper  by  Lamb  (1913)  in  which  the 
disturbance  produced  in  the  flow  of  a  uniform  stream  (of  infinite  depth) 
by  a  submerged  circular  cylinder  is  dealt  with.  The  evaluation  yields. 


actually,  the  disturbance  produced  in  the  flow  of  the  stream  by  the 
dipole  which  generates  the  cylinder  in  an  unbounded  fluid.  The  boundary 
condition  on  the  surface  of  the  cylinder  is  not  satisfied  and  an  anom¬ 
alous  moment  on  the  cylinder  results  (Wehausen  1960,  p.  576;  see  also 
Havelock  1929).  As  Tuck  (1965)  has  shown,  no  closed  body  is  actually 
generated.  Thus,  Lamb's  analysis  provides  only  a  first  step  in  the 
solution  of  the  linearized  problem. 

Let  0 ^  be  the  velocity  potential  of  a  uniform  (unbounded)  stream. 
Let  0 ^  be  the  velocity  potential  of  the  image  of  0 ^  in  the  sub¬ 
merged  body;  that  is,  0 ^  is  such  that  0^  +  0 ^  satisfies  the  boundary 

(s  ) 

condition  on  the  surface  of  the  body.  Let  0  v  '  be  the  velocity  po- 

(s) 

tential  of  the  image  of  0 ^  in  the  free  surface;  that  is,  '  is 

such  that  0 +  0 ^  satisfies  the  free-surface  boundary  condition 

and  the  radiation  condition  at  x  =  -oo.  Let  in  general  0  v  be 

n 

the  velocity  potential  of  the  image  of  0  in  the  free  surface  and 

/g  \ 

0  .  be  the  velocity  potential  of  the  image  of  0  v  J  in  the  submerged 

n+1  n 

body .  The  sequence 

0O  +  01  +  01(S)  +  .  .  .  +  0n  +  0n<s>  +  0n+1  + - 

* 

if  convergent,  yields  the  complete  solution  to  the  linearized  problem. 

This  program  has  been  carried  out  by  Wehausen  (Wehaasen  &  Laitona 
1960,  pp.  574  et  seqq.)  for  the  case  of  the  circular  cylinder.  In 


«•**»?  V3*  * 


this  case,  0  can  be  obtained  f roc  0  by  toeans  of  a  formula  of 

n  n 

'  (s) 

Kochin,  and  0  . ,  from  0  by  using  Milne-Thcmson's  circle  theorem; 

n-*-l  c 

Wehausen  has  shown  the  series  to  be  convergent.  Cylinders  of  other 
shapes  can  probably  be  handled  by  a  combination  of  this  technique  and 
conformal  mapping. 

A  complete  solution  for  the  submerged  circular  cylinder  was  also 

given  by  Havelock  in  1936.  in  previous  papers  (1917,  1927,  1929) 

Havelock  had  applied  to  this  problem  the  method  of  successive  images 

later  used  by  tfehausen  to  give  a  complete  solution,  evaluating  the 

second  set  of  images  within  the  cylinder  and  the  corresponding  image 

in  the  free  surface,  and  thus  carrying  the  computations  two  stages 

further  than  in  Lamb's  solution.  In  his  1936  paper,  however,  he  used 

a  different  approach  and  obtained  the  solution  by  expanding  the  complex 

potential  of  the  system  of  sources  and  sinks  within  the  cylinder  in  a 

-n 

Laurent  series,  Z  A  z  ,  about  the  origin.  The  free-surface  boundary 

n 

condition  (and  the  radiation  condition)  can  then  be  satisfied  by  adding 
a  suitable  expression  in  terms  of  the  coefficients  of  the  Laurent 

series,  and  the  boundary  condition  on  the  surface  of  the  body  yields 
an  infinite  set  of  equations  for  the  coefficients  A^.  The  procedure, 
which  lends  itself  well  to  approximate  computations,  amounts  to  pro¬ 
ducing  the  flow  outside  the  cylinder  due  to  the  system  of  singularities 
inside  by  placing  at  its  center  an  infinite  number  of  multipoles,  mod¬ 
ified  to  satisfy  the  free-surface  boundary  condition  (and  the  radiation 


condition),  and  whose  strength  is  chosen  so  as  to  satisfy  the  boundary 


condition  on  the  surface  of  the  body.  Havelock  found,  for  the  wave 
resistance  of  the  circular  cylinder,  the  expression 


R  «  4s2pU2(k  a)1  e'2k°d 


l-2r  (k  a)2  -  (r,-3r  2+s2)  (k  a)4  +  .  .  .  < 
1  o  2  1  o  r 

j 


where  a  is  the  radius  of  the  cylinder,  d  the  depth  of  submergence  of 
its  axis, 

ko  «  g/u2, 

s  =  2.^e 


a  =  2k  d, 
o 


I 

-  n  -  -  o 

r - —  -r  2 

n  n-rl 

a 


LbiII I  +  IsiilL  + 

n  n- 1 

a  a 


i  -a  ^ 

.  .  t  -  -  e  li(e  ,  ^ 


■i(ea)l  , 


and  li  denotes  the  logarithmic  integral.  The  first  term  in  this  ex¬ 
pression  is  that  obtained  by  Lamb  for  the  wave  resistance  of  the  circu¬ 
lar  cylinder. 

The  approximate  .  to  the  wave  resistance  of  a  submerged  body  which 
is  obtained  by  evaluating  the  flow  about  the  singularity  distribution 
that  produces  the  body  in  an  unbounded  "stream  has  been  calculated  by 
Havelock  for  the  sphere  (1917),  for  prolate  and  oblate  spheroids  moving 
in  the  direction  of  the  axis  of  symmetry  and  at  right  angles  to  it 
(1931a),  and  for  an  ellipso^u  with  unequal  axes  moving  in  the  direction 
of  the  longest  axis  (1931b).  Havelock  evaluated  the  vave  resistance  of 
the  sphere  by  direct  integration  of  the  horizontal  component  of  the 
pressure  on  its  surface;  the  method  of  successive  images  was  applied  to 
compute  the  second  set  of  images  within  the  sphere,  needed  in  order  to 
satisfy  the  boundary  condition  on  its  surface  to  the  required  degree 
of  approximation.  The  computation  of  this  second  set  of  images  can  be 
dispensed  with  if  Lagally's  theorem  is  used  to  obtain  the  W3vi*. 
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resistance  to  the  same  degree  of  approx imat ion.  This  fact  was  noticed 
by  Havelock  in  his  1929  paper  on  the  circular  cylinder,  in  which  a  proof 
of  Lagal'y's  theorem  for  the  steady  two-dimensional  case  in  which  only 
sources,  sinks  and  doublets  are  present  is  given.  The  proof  is  based  on 
the  Blasius  formulas  for  the  force  and  moment  on  a  closed  contour  in 
two-dimensional,  incompressible,  potential  flow;  the  name  of  Lagally  is 
not  associated  with  the  result.  To  compute  the  wave  resistance  for  the 
spheroid  and  the  ellipsoid,  however,  Havelock  made  use  of  a  previous 
result  for  the  wave  resistance  of  an  arbitrary  set  of  doublets  in  a  . 
uniform  stream  (1923:  see  also  1932)  obtained  using  an  artificial  method 
due  to  Lamb  (1926).  The  method  consists  in  making  use  of  Rayleigh's 
artifice  of  including  a  fictitious  dissipative  body  force,  -pp'V, 
proportional  to  the  disturbance  velocity,  calculating  the  energy  dissi¬ 
pation  as 


where  the  integration  is  taken  ever  the  boundary  of  the  fluid,  and 
interpreting  the  finite  limiting  value  of  this  expression,  when  p'  is 
made  to  approach  zero,  as  the  rate  at  which  energy  is  propagated  out¬ 
wards  in  the  wave  motion. 

In  his  1929  paper  on  the  circular  cylinder  Havelock  pointed  out 
also  that  the  second  set  of  images  within  the  cylinder  is  needed  if 
Lagally' s  theorem  is  used  to  verify  (to  the  degree  of  approximation 
afforded  by  the  first  image  within  the  cylinder  and  the  corresponding 
image  in  the  free  surface)  that  the  moment  about  the  center  of  the 
cylinder  vanishes.  Indeed,  application  of  Lagally's  theorem  shows 


immediately  that  a  contribution  of  the  required  order  of  magnitude 


arises  from  the  interaction  of  the  external  uniform  stream  and  the 
second  set  of  images.  For  a  prolate  spheroid  moving  along  its  axis, 
the  contribution  to  the  moment  arising  from  this  interaction  was  calcu¬ 
lated  by  Havelock  (1952).  The  same  task  was  undertaken  by  Pond  (1951) 
for  the  Rankine  ovoid  and  later  (1952)  for  the  more  general  case  of  an 
elongated  body  of  revolution.  Pond  does  not  obtain  the  second  image 
system  exactly.  Instead  he  uses  Munk's  technique  and  obtains,  witn 
some  additional  simplifications,  an  approximate  image  system  in  the 
form  of  a  distribution  of  doublets  along  the  axis  of  the  body  between 
the  limits  cf  the  distribution  that  produces  the  body  in  the  uniform 
unbounded  stream.  Lagally's  theorem  yields  then  immediately  the  re¬ 
quired  moment.  Pond's  result  agrees  with  the  approximation  that 
Havelock,  based  upon  his  analysis  for  tie  prolate  spheroid,  proposes 
for  long  slender  bodies  of  revolution. 

The  case  of  the  fully  submerged  prolate  spheroid,  which  will  be 
dealt  with  in  this  work,  was  treated  analytically  by  Bessho  (1957), 
who  tried  to  satisfy  exactly  the  boundary  condition  on  the  surface  of 
the  spheroid  using  a  distribution  of  sources  on  it.  Bessho  considered 
first  an  ellipsoid  with  three  unequal  axes  and  was  thus  able  to  use  in 
the  solution  of  the  problem  several  results  on  Lame's  functions  given 
by  Hobson  (1931,  Chap.  11).  His  final  expressions,  however,  appear  to 
be  incorrect,  possibly  because  of  typographical  errors,  and  his  numeri¬ 
cal  evaluations  are  rather  inaccurate.  In  the  present  work,  a  more 
direct  attack  using  spherical  coordinates  yields  an  equivalent  set  of 
equations  for  determining  the  source  distribution,  but  with  certain 


significant  corrections.  Furthermore,  the  numerical  evaluations  are 
performed  to  a  high  degree  of  accuracy,  avoiding  the  rough  approxima¬ 
tions  employed  by  Bessho. 
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BASIC  EQUATIONS 

For  convenience  of  reference  the  equations  describing  the 
irrotational  notion,  under  gravity,  of  an  incompressible  fluid  having 
a  free  surface  are  collected  here.  Derivations  of  these  equations 
can  be  found  elsewhere. 

Let  Oxyz  be  a  right-handed  rectangular  coordinate  system  and 
let  the  y  axis  be  directed  vertically  upwards.  Since  the  motion  is 
irrotational,  we  have 

v  «  V<§  (i) 

where  V  is  the  velocity  vector  and  is  the  velocity  potential. 

Since  the  fluid  is  incompressible,  it  follows  from  the  equation  of  con¬ 
tinuity  that  <^>  must  satisfy  Laplace's  equation 

V  2<|  «  o  (2) 

The  Euler  equations  of  motion,  on  the  other  hand,  yield  the  integral 


where  g  is  the  acceleration  of  gravity,  p  is  the  pressure,  p  is 
the  mass  density  of  the  fluid,  and  t  is  time. 

To  these  equations  we  must  add  the  appropriate  boundary  condi¬ 
tions  for  the  problem  on  hand.  Let  S(t)  be  a  boundary  surface  and 
let  it  be  described  by  F(x,y,z,t)=0.  Since  no  transfer  of  matter  takes 


ii rttix  i— niini(i»Twm in— jryjuBBfcfa 
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place  across  S(t),  we  must  have  at  every  point  of  the  surface 


V 

x 


3f 


+  v 
y 


3f 

51 


,  Bf 
+  5E 


o 


(4) 


where  V  ,  V  ,  and  V  are  the  components  of  the  velocity  vector  V. 
x  y  z 

At  a  fixed  boundary,  condition  (4)  reduces  to 

=  0  (5) 


B/dn  denoting  the  derivative  in  the  direction  of  the  outward  normal  to 
the  boundary.  If 


y  =  ^ (x,z,t) 


(6) 


describes  the  free  surface  c£  the  fluid,  condition  (4)  takes  the  form 


V  o)  -  V  +  V  *1  +'>'1  =  0  (7) 

x  (x  y  z  (z  (t 

Let  p^  be  the  constant  pressure  above  the  free  surface.  Then, 
in  addition  to  the  kinematic  condition  (7),  we  must  have  at  every  point 
of  the  free  surface 


p(x,y,z,t)  =  pQ 


(8) 


since  we  neglect  both  surface  tension  and  viscous  effects.  When  use  is 

s 

made  of  Equation  (3),  this  condition  becomes 


(9) 


to  be  satisfied  on  y  =  /,|(x,z,t),  the  additive  constant  having  been 
merged  in  the  value  of  d^/c>t. 

In  the  following,  the  perturbation  produced  by  a  submerged  obstacle 


11 


in  the  flow  of  a  stream  of  constant  velocity  U  will  be  considered.  It 
is  then  convenient  to  write 


i 


=  Ux  +  0 


(10) 


where  the  perturbation  potential  0  satisfies  the  Laplace  equation 


V  0  =  0 


(11) 


and  to  let  u,  v,  and  w  be  the  components  of  the  perturbation 
velocity,  that  is 


d0 

~  dx’ 


(12) 


The  boundary  conditions  become 


U 


0x  c)0 


o 


(13) 


on  the  surface  of  the  body,  and,  since  the  motion  is  steady. 


(U-Hi) 


/yjx  -  v  H-  w  ^  =0 


(14) 


and 


+  l( 


2  2? 
(U+u)  +  v  +  w" 


1  2 
=  2  U 


(15) 


on  the  free  surface  y  =  *|(x,z). 

The  boundary  condition  that  the  velocity  potential  must  satisfy 
on  the  free  surface  can  be  obtained  by  eliminating  ^  between  (14) 
and  (15),  as  was  done  by  Landweber  (1964).  It  is  more  convenient,  how¬ 
ever,  to  proceed  as  follows.  On  the  free  surface 


P(x,y,z)  =  p 


and  therefore 


12 


.  .  ..  .  _  ,  i  ■ 

w 
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that  is,  the  perturbation  vanishes  far  below  the  body.  In  addition  to 
the  boundary  conditions  (13),  (14),  (15),  and  (20)  we  must  also  have 

lim  grad  0=0  (21) 

x — —  -  <=*« 

that  is,  the  motion  must  also  vanish  far  ahead  of  the  body.  This  so- 
called  "radiation  condition"  need  not  be  imposed  if  Rayleigh's  artifice 
of  introducing  a  fictitious  dissipative  body  force  proportional  to  the 
disturbance  velocity  is  used  (Lamb  1932,  p.  399;  see  also  Wehausen  and 
Laitone  1960,  p,  479)  or  if  an  initial  value  problem  can  be  formulated 
(for  which  the  boundary  condition  at  infinity  is  simply  that  the  motion 
vanishes  everywhere)  whose  solution  yields  the  steady-state  solution 
sought  as  t— —  (Stoker  1957,  pp.  174-181;  see  also  Wehausen  and 

Laitone  1960,  p.  472,  for  additional  references). 

Although  Laplace's  equation  (11)  and  the  boundary  condition  (13)  on 
the  boundary  of  the  body  are  linear,  the  free  surface  boundary  condi¬ 
tions  (14)  and  (15)  are  nonlinear.  This  nonlinearity  precludes  the  use 
% 

of  techniques  of  solution  based  on  the  principle  of  superposition. 

Thus,  the  method  of  separation  of  variables  and  expansion  in  eigen¬ 
functions  cannot  be  used,  and  neither  can  the  method  of  Green's  func¬ 
tions  or  singularity  distributions.  Moreover,  the  free-surface  boundary 
conditions  are  rather  inconvenient  because  they  are  to  be  applied  on  a 
surface  of  unknown  location. 

Under  these  conditions,  special  problems  in  jhe  theory  of  surface 
waves  have  been  treated  by  using  approximation  techniques  (actually 
perturbation  procedures)  of  which  a  fairly  complete  account  has  been 
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given  by  Wehausen  (Wehausen  and  Laitone  1960).  Since  only  the  first 
order  approximation  will  be  treated  here,  the  perturbation  analysis 
applicable  in  our  case  is  r^t  included  and  the  reader  is  referred 
to  Wehausen1 s  treatise  for  details  of  the  method.  The  first  order 
approximation  consists  of  neglecting  all  second  order  terms  in  the 
free-surface  boundary  conditions  and  applying  them  on  the  plane  y  =  0 
instead  of  on  the  free  surface  y  =  '’'[(x.z).  Equations  (14),  (15),  and 
(19)  become,  respectively. 


U^x  -  v  =  0  (22) 

g  +  Uu  =  0  (23) 


and 


(24) 


to  be  satisfied  on  the  plane  y  =  0.  Equation  (24)  can  be  obtained 
directly  from  (22)  and  (23)  by  eliminating  between  these  two  equa¬ 

tions  . 


It  is  interesting  to  note  the  boundary  condition  that  the  second- 
order  contribution  to  the  perturbation  potential  must  satisfy.  This  is 


(25) 


to  be  satisfied  on  the  plane  y  =  G,  Once  the  first-order  approximation 
0^^  is  known,  therefore,  the  second-order  contribution  -:an  ae  obtained 
by  solving  another  linear  problem  with  a  nonhcraogeneous  boundary  con¬ 
dition  on  the  plane  y  =  0. 


- nr—rrumwwM*, 
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ANALYTICAL  TREATMENT 


The  solution  will  be  sought  in  the  form  of  a  distribution  of 
sources  on  the  surface  S  of  the  spheroid.  The  determination  of  the 
surface  density  a  of  this  distribution  is  the  object  of  the  follow¬ 
ing  analysis. 

Let  the  x  axis  coincide  with  the  major  axis  of  the  spheroid 
situated  at  a  depth  d  below  the  free  surface.  Choose  the  y  axis 
vertical  and  upwards  and  let  the  origin  be  at  the  center  of  the  sphe¬ 
roid.  Let  the  stream,  of  constant  velocity  U,  flow  in  the  positive 
x  direction.  The  velocity  potential  corresponding  to  the  three- 
dimensional  motion  past  a  unit  source  at  )  is  given  by 

^s  =  ^  "  R  +  G(x,y,z;^ ’7  (26) 


where 


6-(»/3,a;X7,^)=  ^  I  ™ I  aZ 


2*r  ,oo 


Uo  sect 


zu  -k(d-y) 


ko  sectt 


dk  dt 


4 


£ 

Z 


"Z  /<*,  l  I  Sec  t 


t  (C^~7JJ  /c0  %ec  +6  ks  t -b  tj 

e  e  ^ 


-H 

z 


dt 


(*7) 


y  <  2d  - 
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*  2 

Here  Re  denotes  the  real  part  of  a  complex  number,  =  g/U  ,  g  is 

the  acceleration  of  gravity,  and 


R2  =  (x-  £)2  +  (y-''])2  +  (z-)  )2 

The  harmonic  function  G(x,y ,z ;  5 ,  "|  , ^  )  is  sometimes  called  the 
Havelock  potential  since  Havelock  (1932)  was  the  first  one  to  express 
it  as  a  double  integral  in  terms  of  Rayleigh's  fictitious  viscosity. 
We  have  then  for  the  velocity  potential  c]>  corresponding  to  the  dis¬ 
tribution  of  sources  on  the  surface  S  of  the  spheroid  the  express icr 

's  's 

y  <  2d  -  b 


where  2b  is  the  length  of  the  minor  axis  of  the  spheroid  and  the 
integrations  over  the  surface  S  of  the  spheroid  are  performed  ir, 
the  variables  \  ,  and  ^  . 

Since  (26)  satisfies  the  free-surface  boundary  condition,  the 
boundary  condition  at  y  =  -  00.  and  the  radiation  condition  at 
x  =  _  00 ,  so  does  the  velocity  potential  given  by  (28).  The 

unknown  surface  density  a  must  therefore  be  determined  so  that  the 
boundary  condition  the.  surface  of  che  body. 


The  symbol  Re  could  be  left  out  since  the  imaginary  parts  of  the 
two  expressions  within  brackets  are  identically  zero. 
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is  satisfied.  VJhen  <j>  ,  as  given  by  (28),  is  substituted  into  this 
boundary  condition,  a  Fredholm  integral  equation  of  the  second  kind  is 
obtained  for  c, 

Z"-cU7,*)+j  nVAVl^j-^  +■  =  <29> 

The  solution  of  this  integral  equation  will  be  sought  by  expanding  the 
fimctions  involved  in  series  of  spheroidal  harmonics,  thus  taking  ad¬ 
vantage  of  the  particular  shape  of  the  body. 

On  changing  the  order  in  which  tha  integrations  are  performed,  we 
can  write  in  the  equivalent  form 


Ux-|  ~ 

izj  +  t  Ssu.-tj 


LOT'o 


21 T  ,  °°  t  ; 

,(  ;zh‘L 


hi 

Si 


k’~  c* 


j  e  ^(V/ZVH3 

S  J 


3L 

'2 


/2  ~zk^i  s-ec  t  kcs&cr\  ^ (x.<bcf+£sii<t}l 

2k*L  j  s#c  t  -e.  €  L 

iff 

2. 


l 


•SacT  t  J"  ^  -  t  (  ^  <kS  ^  Sil.f'jj 


k.  -  * 


(30) 


s 

We  will  denote  the  length  of  the  major  axis  of  the  spheroid  by  2a  and 
its  focal  distance  by  c.  Let  ^  ,  &,  and  be  the  orthogonal  con- 
focal  coordinate  system  defined  by 

x  =  c  cosh  ^  cos  G  , 

y  =  c  sinh  ^  sin  ©  cos  , 

s  inh  'J'j  s  in  9  s  in  ^  . 


z  =  c 


(31) 
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Th  surfaces  =  constant  are  confocal  prolate  ellipsoids  of 
revolution,  with  common  foci  at  the  points  (+c,0,0).  The  surfaces 
0  =  constant  are  confocal  hyperboloids  of  revolution,  with  the  same 
foci.  Let  q  be  the  value  of  the  confocal  coordinate  that 

corresponds  to  our  spheroid.  We  have 


a  =  c  cosh  W  , 


b  =  c  sinh 


Let  f(.9  be  the  value  of  the  Newtonian  potential 

'  s  ‘ 

on  the  surface  of  the  spheroid,  and  let  it  be  represented  by  a  series 

* 

of  surface  spherical  harmonics.. 


Hobson's  notation,  which  is  also  the  one  chosen  by  the  National  Bureau 
of  Standards  in  its  Handbook  of  Mathematical  Functions,  is  used  in  this 
work.  We  have 

m 

m  9  2  dmP  (z) 

P  (z)  =  (z  -1)  - —  , 

a  .  ,  \  m 

(dz) 

m 

o  2  dmP  (x) 

P„(x)  (-1)  (1-x  )  -  » 

n  n  \m 

(dx) 

m 

„  2  dmQ  (z) 

Q  (z)  =  (z  -1)  - —  , 

^  f  A  \ ^ 

(dz) 

m 

„  2  dmQ  (x) 

Q  (x)  =  (-1)  (1-x  ) - —  , 

n  ,  ,  Nm 

(dx) 

where  z  is  any  complex  number  not  real  and  between  -1  and  1,  x  is 
a  real  number  in  the  interval  -1  <  x  <  1,  n  and  m  are  positive  inte¬ 
gers  or  zero,  m  <  n.  The  only  points  of  discontinuity  of  the  functions 
Pn(z)  (m  i-  0)  and  Qn(z)  are  on  the  straight  line  (-1,+1). 
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l 


M:Q  A»jr  O 


We  assume  this  scries  to  be  absolutely  and  uniformly  convergent  over 
the  surface  S.  We  can  then  write  (Hobson  1931) 

o\ 

\ —  r —  /*w 

/  =  -/  ”tS  =  )  )  <T  P>se)^  04) 

'  *-  L-  Q"VM) 

■n-O  w=c>  vi  V  io; 

this  being  the  potential  function  for  the  space  external  to  the  spheroid 
which  converges  to  f  ( O ,  ip" )  on  the  boundary.  For  the  potential  func¬ 
tion  for  the  internal  space  which  converges  to  f over  the 
boundary  we  have  the  corresponding  series  expansion 

o O  VI 

y  y  (35) 


The  surface  density  a  can  be  expressed  in  terms  of  the  coefficients 
a^  of  these  expansions.  We  obtain  (see  Appendix  A) 


Cv  (V) 


=  u 


AtcO  Avf 


“  <W1  5 

A  -U1  - _ (36) 

ATT 


where,  for  convenience,  we  have  put 


a  1 

A  —  a 

»nt  * 
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<-  v  *  rfs*  ' 


arsaMB’ 


Substituting  this  expression  for  CT  in  the  surface  integrals  contained 
in  the  third  and  fourth  terms  on  the  right  side  of  (30),  we  find  for 
the  first  of  these 

/  «  <r(^)As  =^r)  )  H) 


i  A  )  +■ - ! _ 


z  Z_  L 

*\  :r  O  'yv\  tc  O 

\pL  J  IA) 
jzA,  Vr^'’ 


(38) 


where 


A  =  k  c  cos  t. 


since  (see  Appendix  C,  Equation  22  with  i  replaced  by  -i) 

£  ‘  \Co--S) 


d  - 


Cs  r.tU^(feA2yo-4/&)^ 


TT/  ^1/  Lcclf’ -  <-  Ai©(crSr+4l^tM  0S«m(P 

0  ’0 


Su,  &  cl&d 


«i+-iwi  m+wI  am  /  \ 

=  2!r(-')  ^  'IMlo)' (**t+fc-.t-)  +  — — — - 

s  pccr+^t-j  j 


;X  J  (A)  (39) 

■!2A  «Hf 


(We  assume  that  the  integration  can  be  carried  out  term  by  term.)  The 

surface  integral  contained  in  the  fourth  terra  on  the  right  side  of  (30) 

2 

is  given  by  (38)  with  k  replaced  by  k^sec  t. 

Finally,  using  (38)  and  the  corresponding  expression  for  the  sur¬ 
face  integral  in  the  fourth  term,  and  expanding  the  functions 

2-f . :)] 


ky+ik(x  cos  t  +  s  sir,  t) ,  eko«c't  [y+i(x  cos  t  +  z  sin  t: 


in  series  of  spheroidal  harmonics,  we  can  express  the  velocity  poten¬ 
tial  <P  in  the  form 


>a&a&*»£se 


1 
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4-^e 


TT 
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(40) 


where  the  coefficients  ^  are  given  by  (24)  of  Appendix  C,  with 

C  =  C  ,  the  coefficients  B  have  been  omitted  because  the  po- 
n,o  n  n,m  r 


tential  is  an  even  function  of  z,  the  coefficients  C 


are  ob- 


n,m 


tained  from  the  coefficients  C  by  putting  k  =  k  sec^t,  and, 

n,m  J  r  6  o  *  * 

for  brevity,  the  functions 


=  ±|(sech-k-^  +■ 


X  l/tc4>sfc") 

z 


(41) 


and 


T  m(t)  =  Tm(k  sec2t,t) 
o  n  o 

n 

have  been  introduced.  Again  assuming  that  integration  of  the  series 

term  by  term  is  possible,  and  replacing  the  coefficients  C  and 

n,m 

C  by  their  expressions  in  terms  of  the  functions  Tm  and  T  m, 

o  no’ 

n,m  n 


J 

C  -  , 
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we  can  write  (40)  in  the  equivalent  form 
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^  i-  ^  ^i' 

where  the  factor  1/2  is  to  be  used  for  m  -  0.  Substitution  of  this 
expression  for  the  velocity  potential  into  the  boundary  condition  on 


the  surface  of  the  spheroid, 

5n  - 


0  or 
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******* 


whera  we  have  now  assumed  that  the  series  can  be  differentiated  term 
by  term;  the  dot  above  the  Legendre  functions  P  and  Q  indicates 
differentiation  with  respect  to  the  argument  cosh  "j  .  Equation  (43) 


is  equivalent  to  the  infinite  set  of  equations 
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CALCULATION  OF  THE  WAVE  RESISTANCE 


Lagalxy's  theorem  (Lagally  1922)  yields  for  the  norizontal  component 
of  the  force  on  the  spheroid  the  expression 


—  tS 


where  is  the  Havelock  potential  for  the  source  distribution  on  the 

surface  of  the  spheroid,  given  by  the  third  and  fourth  terms  on  the 
right  side  of  Equation  (30).  We  have 
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where  both  surface  integrations  are  performed  over  the  surface  of  the 
spheroid,  in  the  variables  x,  y,  and  z  when  we  denote  it  by  S,  and 
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in  the  variables  ^  ,  and  ^  when  we  denote  it  by  S’.  The 

contribution  to  the  wave  resistance  from  the  first  term  on  the  right 
side  of  (46)  is  zero.  To  see  this  we  need  cnly  write 
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and  observe  that  if  we  interchange  x,  y,  z,  with  ,  ^/  >  \  >  the 

integrand  changes  sign  while  its  absolute  value  remains  unchanged.  To 

evaluate  the  contribution  to  the  wave  resistance  from  the  second  term 

on  the  right  side  of  (46),  we  make  use  of  (38)  with  k  replaced  by 
2 

k  sec  t.  We  obtain 
o 


^4* 


cr  ds  -  & 

0* 


~zkol 


JL 


"Sec.  t  -€ 


C  k0  t  [krd-  j 


's  ^  L_  ~ZL 

?  | 

oo  OO  /VI  ^  - 

\  \  'OI  +  '  /v»'i  _  »V»  1  f  \  \  *'*+■•  ,1  y-rr-  —2™ 

Cl  "  -  - 


fflrO 


11  =  0  Mt'zO 


The  wave  resistance  is  then  given  by 

t 


0^2  /T)  OO  rrt 


,  I  tm  rtn  I  >m  »» 


t-0  A*, 


(47) 


if =0  fkt-O  12-0  >rn^O 

(n+m+n '+m' )  even 
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NUMERICAL  ~ VA:  .UAT IONS 


All  numerical  evaluations  were  performed  with  an  IBM  360/65  com¬ 
puter-  The  corresponding  programs  can  be  found  in  Appendix  E  and  will 
be  discussed  in  this  section.  Subprograms  which  are  common  to  various 
programs  appear  only  once  in  the  Append ix. 

In  order  to  evaluate  the  single  integrals  in  the  coefficients  of 
the  infinite  system  of  equations  (44),  we  introduce  the  change  of 
variable 


tan  t  =  u. 


Let 


-L, 


s ecztr 

Te„(t)To  l t)dt 

’  Q  CA-X—»  l.  ^ 


A  -zk0A 

=  -8TT  ^ - 

I  o  4oc>  t 


We  have 


-2.!^  ( 


T,=  « 


00  -ZLduz 


\J  M-  u: 


o  (yjTt;  luHu 

m  7  ai/ 


where 


Ly)= 1 1  ,,X  >4 


(48) 


(49) 


[  (^LH^+upJ  «+" i 

Since  the  integrand  in  (49)  is  an  even  function  of  u,  we  can  write 


(50) 


whers  we  have  further  changed  the  variable  of  integration  to 


x  =  \|  2k  d  u 

1  o 


The  integral  in  (50)  can  be  evaluated  by  Hermite -Gauss  quadrature.  A 
twenty-point  formula  was  used,  for  which  the  zeros  and  weighting  factors 
can  be  found  in  Abramowitz  and  Stegun  (1964).  Actually,  since  the 
integrand  is  an  even  function  of  x,  only  ten  ordinates  are  evaluated 
in  the  corresponding  c  mputer  program,  which  yields  the  value  of  1^ 
and  includes  a  function  subprogram  (BESSJ)  to  evaluate  the  Bessel 
functions  which  appear  in  the  integrand. 

The  evaluation  of  the  Bessel  functions  J  of  order  n  +  1/2  was 
carried  out  using  the  series  expression 


T  =  - )  | - Si -  4.  - “ -  (5] 

r("+-4)  {  z(2"+i!}  2.4'(2^fS)( 2^+5)  j 

for  values  of  the  argument  less  than  five,  and  the  recurrence  relation 


t,  e>  =  r  x  t 


for  values  greater  than  five.  In  this  manner,  the  large  round-off 
errors  associated  with  the  use  of  (51)  for  large  values  of  the  argument 
(especially  for  small  n) ,  and  with  the  use  of  (52)  for  small  values 
of  the  argument,  are  avoided.  The  cut-off  value  of  five  is  a  rough 
estimate  of  the  value  necessary  to  minimize  these  errors.  The  sub¬ 
program,  which  was  written  in  doubl-  precision  further  to  reduce  round¬ 
off  errors  (only  six,  or  at  most  seven,  significant  digits  are  supplied 
by  the  IBM  360/65  computer  in  single  precision),  was  checked  against 
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tables  included  in  Watson's  treatise  (Watson  1944,  pp.  740-743).  The 
results  were  found  to  be  accurate  to  the  number  of  significant  figures 
given  in  these  tables  (six  or  less)  for  values  of  n  up  to  18;  most 
of  them  are  probably  accurate  to  several  more  significant  figures.  The 
round-off  errors  associated  with  (52),  however,  increase  with  increasing 
n,  and  this  high  accuracy  should  not  be  expected  to  obtain  in  the  eval¬ 
uation  of  Bessel  functions  of  much  larger  order. 

In  view  of  the  difficulty  of  estimating  the  error  in  the  Hermite- 
Gauss  quadrature  formula,  the  calculations  were  spot-checked  using 
Simpson's  rule,  which  requires  longer  computing  times,  but  allows  for 
a  ready  estimate  of  the  error  at  each  stage  of  apr roximation .  The  sub¬ 
routine  SMPSN  used  in  the  corresponding  program  is  part  of  the  IBM 
System/360  Scientific  Subroutine  Package,  Version  II,  and  is  not  in¬ 
cluded  in  Appendix  E.  In  this  subroutine,  Simpson's  rule  is  used  with 
interval  halving  until  the  difference  between  successive  values  of  the 
integral  is  less  than  a  given  tolerance.  The  Hermita-Gauss  quadrature 
formula  was  found  to  provide  at  least  five  significant  figures. 

To  evaluate  the  double  integrals  in  the  coefficients  of  (44),  we 
make  use  again  of  the  change  of  variable 

tan  t  =  u 


and  we  put  further 


k  =  k  K 
o 


Let 


.ST 


I*-*'™ 
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We  have  then 
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-h'  f”2 i d kc(u  <54) 

For  nuir.erical  evaluation,  it  is  better  to  write  (54)  in  the  f  rm 


X-7  =r  A'0C 
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(55) 


where  the  function 


>  i 

S*  4.  1 

XT- 

2-i  Pr+!)  [ 

2fZ^+=>)  J 

continuous  at 

z  =  0.  A  simple 

modification  of  the  function 

(56) 


program  used  to  evaluate  the  Bessel  functions  yields  the  function  sub¬ 
program  (BESSJM)  used  to  evaluate  (56). 

The  evaluation  of  the  Cauchy  principal-value  integral  in  (55)  was 


«* ******  ■ 


*W&W  a 
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carried  out  using  Simpson's  rule,  modified  according  to  a  technique 
suggested  by  L.  T.andweber  (Kobus  1967,  Appendix  1)  to  take  into  account 
the  singular  point  in  the  integrand.  It  can  be  shown  that  any  composite 
integration  rule,  obtained  by  dividing  the  interval  of  integration 
into  subranges  and  applying  in  each  subrange  any  of  the  Newton-Cotes 
integration  formulae,  can  be  used  to  evaluate  a  Cauchy  principal-value 
integral, 

a  <  c  <  b 
a 

as  long  as  the  singular  point  c  is  either  an  end  point  (excluding  a 
and  b  )  or  the  middle  point  of  a  subrange,  and  the  value  of  the  inte¬ 
grand  at  c  is  taken  equal  to  the  derivative  f'(c).  In  the  function 

subprogram  that  evaluates  the  function  G(u)  defined  in  (55),  there- 

2 

fore,  the  derivative  at  the  singular  point  k  =  1+u  is  obtained  first 
% 

using  a  five-point  differentiation  formula.  The  two  function  subpro¬ 
grams  FM  and  F  evaluate  the  integrand  of  the  principal-value  inte¬ 
gral  and,  at  the  singular  point,  replace  the  integrand  with  the  corres-. 
ponding  derivative.  The  calculations  were  carried  out  to  an  accuracy 
of  three  significant  figures. 

The  solution  of  the  infinite  system  of  equations  (44)  was  obtained 
by  using  the  so-called  method  of  reduction  (Kantorovich  and  Krylov  1958, 
pp.  25-26,  30-31).  In  th,is  method,  the  solution  is  found  by  solving  a 
sequence  of  finite  systems,  each  of  which  is  obtained  from  the  infinite 
system  by  discarding  all  equations  and  unknowns  beyond  a  certain  number 
N.  The  solutions  of  this  sequence  of  systems,  under  certain  conditions, 
converge  to  the  (principal)  solution  of  the  infinite  system.  Although 


I 


it  is  difficult  in  our  case  to  show  that  the  conditions  of  the  theorem 
are  satisfied,  it  is  reasonable  to  assume,  on  account  of  the  significance 
of  (44),  that  the  method  does  yield  the  solution  sought.  A  numerical 
verification  of  this  assumption  is  obtained  if,  when  the  number  N  of 
equations  kept  is  increased,  the  coefficients  A™  approach  limiting 
values . 

Such  a  numerical  verification  has  been  included  in  the  computer 

program  for  the  solution  of  (44).  The  number  of  equations  has  been 

taken  equal  to  2,  5,  9,  14,  and  20,  equivalent  to  keeping,  in  expansions 

(34)  and  (36),  only  those  terms  containing  Legendre  functions  of  degrees 

up  to  1,  2,  3,  4,  and  5,  respectively.  The  program  contains  a  sub- 

in 1 

routine  subprogram  for  computing  the  coefficient  matrix  B  ,  from  the 

n  n 

•  m  m f 

symetric  matrix  formed  by  the  double  and  single  integrals  I  ,  .  No 

n  n 

table  was  found  to  give  the  values  of  the  Legendre  functions  and  their 
derivatives  in  the  interval  of  interest  to  our  work,  and  their  evalua¬ 
tion  is  carried  out  in  a  second  subprogram.  The  derivatives  occur  in 
in  m1 

the  coefficients  B  ,  ;  the  Legendre  functions  of  the  first  kind  in 
n  n 

expansion  (36).  No  use  is  made  in  this  second  subprogram  of  the  re¬ 
currence  relations  for  the  Legendre  functions,  on  which,  for  example, 
the  tables  published  by  the  National  Bureau  of  Standards  (Lowan  1945) 
are  based;  severe  round-off  errors  were  found  to  be  associated  with  the 

use  of  these  relations,  in  particular  near  x  =  1.'  Instead,  since  only 

* 

the  functions  of  deg  >  es  up  to  five  were  needed,  the  explicit  expres¬ 
sions  of  these  functions  were  derived  and  used  in  the  computations. 

These  expressions  were  checked,  using  the  National  Bureau  of  Standard 
tables,  for  several  values  of  the  argument;  all  significant  figures 


'wrin--‘T~mirr  rrrrnrrwrTrpii  ini 


given  in  the  tables  were  found  to  be  correct.  The  subroutine  GELG 
u&ed  in  the  program  for  the  solution  of  (44)  is  part  of  the  IBM  System/ 
360  Scientific  Subroutine  Package,  Version  II,  and  is  not  included  in 
Appendix  E.  In  this  subroutine,  the  solution  of  a  system  of  general 
simultaneous  linear  equations  is  obtained  by  means  of  Gauss  elimination' 
with  complete  pivoting. 

Once  the  coefficients  A™  are  known,  we  can  evaluate  the  density 

n  7 

of  the  source  distribution  making  use  of  (36).  A  program  was  prepared  • 
to  obtain  plots  of  this  density  along  meridian  curves  of  the  spheroid. 
The  subroutine  PNMXG1  used  in  this  program,  which  evaluates  the 
Legendre  associated  functions  of  the  first  kind  for  values  of  the  argu¬ 
ment  greater  than  one,  is  simply  part  of  the  subroutine  LEGF;  the 
subroutine  PNMXL1,  which  computes  these  functions  for  values  of  the 
argument  less  than  one,  was  obtained  by  introducing  the  necessary 
changes  iu  sign  in  the  expressions  in  the  subroutine  PNMXG1.  Neither 
subroutine  is  included  in  Appendix  E. 

The  evaluation  of  the  wave  resistance  is  carried  out  in  the  V.st 
program  in  Appendix  E.  A  simple  modification  of  the  program  for  the 
evaluation  of  the  single  integral  in  the  coefficients  of  (44)  yields 
the  program  for  evaluating  the  single  integrals  in  the  expression  (47) 
for  the  wave  resistance:  the  modified  program  is  not  included. 
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RESULTS  AND  DISCUSSION 

The  double  and  single  integrals  mim,  in  the  coefficients  of 

n  n 

system  (44)  depend  on  the  Froude  number  U/ \r2gc  =  l/N^k^c,  and  on 

the  relative  depth  of  submergence  d/c,  and  are  independent  of  the 

eccentricity  c/a  of  the  spheroid,  which  occurs  in  the  derivatives  of 

the  Legendre  functions.  This  property  of  the  coefficients  of  system  (44) 

is  of  interest  because  most  of  the  computer  time  required  in  order  to 

evaluate  the  coefficients  Am  is  spent  in  the  evaluation  of  the  double 

n 

integrals;  once  these  integrals  have  been  computed,  for  given  values  of 
the  Froude  number  and  the  relative  depth  of  submergence,  the  coefficients 
Am  can  be  obtained  for  several  values  of  the  eccentricity  without 
practically  increasing  the  total  computing  time.  The  numerical  evalua¬ 
tions  reported  here  were  carried  out  for  a  Froude  number  U/ \  2gc  =  0.4, 

\ 

a  relative  depth  of  submergence  d/c  =  0.5,  and  reciprocal  eccentrici¬ 
ties  a/c  =  1.01,  1.02,  1.04,  1.06,  1.08,  and  1.10,  corresponding  to 
slenderness  ratios  a/b  =  7.12,  5.07,  3.64,  3.02,  2.65,  and  2.40,  re¬ 
spectively.  (For  the  relative  depth  of  submergence  chosen,  the  spheroid 
pierces  the  free  surface  for  a/c  =  \|l .25 . )  The  practical  interest  of 
the  smaller  slenderness  ratios  is  probably  limited;  they  are  included 
here  for  purposes  of  comparison  and  discussion  of  results. 

The  first  approximation  to  the  solution  of  system  (44)  which  con¬ 
sists  of  neglecting  entirely  the  infinite  series  on  the  right  side,  thus 
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keeping  only  the  first  term  of  expansion  (34)  with 


U  c 


a/c  _  _l  Q-/C-H 


(57) 


corresponds  to  the  source  distribution  that  produces  the  spheroid  in  an 
unbounded  uniform  stream  without  a  free  surface.  Indeed,  the  potential 
for  the  motion  of  a  prolate  spheroid  parallel  to  its  axis  (in  the  nega¬ 
tive  x  direction)  in  an  infinite  mass  of  liquid  is  given  by 


j>  =•  U  c.  — j -  fj5  (ks  o)  0.^-L  'r|  J 


(58) 


or 


j>=  ^  u - x. 

X-'  z 

with  ^  Q  =  cosh^  =  a/c,  =  coshi|  ,  p  =  cos  0.  For  this  first 
approximation,  the  wave  resistance  is  given  by 


7c  -  Tq  c°(A ,°)  ^Tf  j  *eite 


-zkact  5<?<r  r 


4  ^ 


(59) 


Equation  (59)  was  obtained  by  Havelock  (1931a)  using  the  axial  source 
distribution  corresponding  to  the  motion  of  the  spheroid  in  an  infinite 
mass  of  liquid  and  applying  Lagally's  theorem  to  evaluate  the  wave  re¬ 
sistance  (without,  however,  mentioning  the  name  of  Legally  in  this 
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connection).  Of  course,  neither  the  axial  source  distribution  nor  the 
source  distribution  on  the  surface  of  the  spheroid  actually  produce  a 
spheroid  in  the  presence  of  a  free  surface.  They  are  two  different 
first  approximations.  It  is  interesting  to  note,  however,  th.at  both 
first  approximations  yield  the  same  value  as  a  first  approximation  for 
the  wave  resistance.  Havelock's  approximation  has  been  included  in 
Tables  1  through  5  for  purposes  of  comparison. 

The  values  of  the  coefficients  A™,  obtained  keeping  1,  2,  5,  9, 

14,  and  20  equations  (and  an  equal  number  of  coefficients)  in  (44), 
are  presented  in  Tables  1  through  4,  for  four  values  of  the  ratio  a/c. 
Table  5  contains  the  corresponding  successive  approximations  to  the  wave 
resistance  (the  two  additional  a/c  ratios  investigated  are  also  in¬ 
cluded  in  this  table).  Three  significant  figures  are  given  for  all 
elements  in  the  tables,  since  the  double  integrals  were  computed  to  this 
accuracy.  When  20  equations  are  kept,  92  double  integrals  must  be  evalu- 

ated  (the  integrals  I  .  satisfy  the  symmetry  relations  I  ,  -  ,1 

nn  J  J  nn  nn 

and  mlm,  =  m  I  m. ,  and  this  reduces  the  number  of  double  integrals  to 
n  n  n  n 

be  evaluated  from  202  to  92);  the  required  computing  time  is  a  little 
less  than  two  hours.  The  corresponding  computing  time  for  the  evalua¬ 
tion  of  the  single  integrals  is  less  than  half  a  minute;  for  the  solution 
of  (44),  for  the  six  ratios  a/c  investigated,  including  the  evaluation 
of  the  derivatives  of  the  Legendre  functions  which  appear  in  the  coeffi¬ 
cients  of  the  system,  less  than  one  minute.  Also  less  than  one  minute 
is  the  computing  time  required  for  the  evaluation  of  the  wave  resistance, 
using  (47);  only  single  integrals  appear  in  this  expression.  A  reduction 
in  the  time  required  for  the  evaluation  of  each  double  integral,  at 
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Table  5.  Wave  resistance  R/jtpgc 
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present  about  75  seconds,  would  therefore  reduce  practically  in  propor¬ 
tional  form  the  total  computing  time.  This  reduction  is  believed  possi¬ 
ble  and  is  part  of  a  future  program  of  work. 

As  should  be  expected,  the  closer  the  spheroid  is  to  the  free 
surface,  the  larger  the  number  of  coefficients  that  must  be  kept  in 
order  to  obtain  the  limiting  values  to  a  given  accuracy.  Most  likely, 

m  I  I 

not  all  products  A  ,  in  (44),  corresponding  to  a  given  number  of 

equations,  need  be  kept  in  order  to  obtain  the  solution  to  that  same 
accuracy.  Moreover,  not  all  coefficients  A™,  up  to  a  certain  one  in 
the  order  in  which  they  appear  in  expansion  (34),  need  be  computed  in 
order  to  obtain,  for  example,  the  wave  resistance,  also  to  that  same 
accuracy,  since  not  all  terms  in  the  summation  in  (47)  contribute 
significantly  to  the  final  value.  It  is  difficult,  however,  to  esti¬ 
mate  beforehand  the  error  associated  with  such  approximations,  a  point 
of  practical  interest  since  a  significant  reduction  in  computing  time 
could  result  thereof;  no  attempt  has  been  made  in  this  work  to  obtain 
such  estimates.  It  is  interesting  to  note  that  when  only  A^  is  kept 
(one  equation),  despite  the  fact  that  the  resulting  approximation  for 
A^  is  close  to  the  limit  value,  the  correction  for  the  wave  resistance 
is  rather  small,  since  it  turns  out  that  some  of  the  coefficients  left 
out  contribute  rather  significantly  to  its  final  value.  Moreover,  keep¬ 
ing  A^  in  addition  to  A^,  does  not  improve  the  approximation  to  the 
wave  resistance,  whicii  remains  unchanged  to  the  three  significant  fig¬ 
ures  given  in  Table  5. 

A  set  of  equations  essentially  equivalent  to  (44)  for  determining 
the  density  ^f  the  source  distribution  on  the  surface  of  the  spheroid 


I 


44 


was  obtained  by  Bessho  (1957),  using  an  entirely  independent  derivation. 
Bessho  considered  first  an  ellipsoid  with  three  unequal  axes,  and  was  thus 
able  to  use  in  his  solution  several  results  on  Lame'f-  .ictions  given  by 
Hobson  (1931,  Chap.  11).  The  coefficients  of  his  set  of  equations,  how¬ 
ever,  appear  to  be  incorrect,  possibly  because  of  typographical  errors, 
and  his  numerical  evaluations  are  rather  inaccurate.  His  approximation 
to  the  solution  of  the  infinite  set  of  equations,  which  contains  only  six* 
of  the  coefficients  of  the  system,  all  belonging  to  the  first  five  equa¬ 
tions,  introduces,  in  the  light  of  the  present  numerical  results,  a  sig¬ 
nificant  error.  Moreover,  rather  than  using  the  singular  solution  for  a 
source  in  the  form  of  a- double  integral  and  a  single  integral  as  done  in 
the  present  work,  the  coefficients  of  the  infinite  system  of  equations 
are  expressed  in  terms  of  Rayleigh's  "fictitious  viscosity"  p,  and  are 
not  therefore  in  a  form  suitable  for  numerical  evaluation;  the  calcula¬ 
tions  are  performed  using  asymptotic  expansions,  with  which  large  errors 

C 

are  associated  unless  the  Froude  number  is  small  enough.  For  a  Froude 

number  U/  ^ 2gc  =  0.395,  the  relative  depth  of  submergence  used  in  the 

present  calculations,  d/c  =  0.5,  and  a  ratio  a/c  =  ^17/ 4  ^1.03,  Bessho 
2  2 

obtains  R/4npU  c  =  7.165x10  ,  a  1467,  increase  over  the  value  corres¬ 

ponding  to  Havelock's  approximation,  2.91x10  ^ .  For  Froude  number  0.4, 
Table  5  shows  a  34%  increase  over  the  value  corresponding  to  Havelock's 
approximation  for  a/c  -  1.02,  and  a  907,  increase  for  a/c  =  1.04. 

The  relative  density  a/U  of  the  source  distribution  on  the  sur¬ 
face  of  the  spheroid,  corresponding  to  the  successive  approximations 
obtained  keeping  increasing  numbers  of  equations  in  (44),  for  a/c  =  1.02, 
is  presented  in  Figs.  1,  2,  and  3.  These  figures  show  that,  although 
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Figure  2. 


Density  of  source  distribution  along  a  vertical  meridian 
plane  near  rear  of  spheroid  (a/c  =  1.02) 


LATIVE  DENSITY  <r/U 


ANGLE  6  (DEGREES) 

si  try  of  source  distribution  along  meridian  planes  of  the  spheroid  (a/c  -  1.02) 


for  an  integrated  quantity  like  tha  wave  resistance,  twenty  equations 
are  enough  to  obtain  the  solution,  the  convergence  is  slower  for  the 
surface  density.  Since  the  limiting  values  of  most  of  the  coefficients 
A™  in  Table  2  have  been  practically  attained  (with  the  probable  excep¬ 
tion  of  those  in  the  last  row),  this  implies  that  some  of  the  coeffi¬ 
cients  left  out  have  a  significant  effect  on  the  value  of  the  density 
of  the  source  distribution.  It  is  interesting  to  note  that  the  surface  - 

t 

density  oscillates  about  the  final  value;  the  approximation  obtained 
keeping  five  equations  in  (44)  is  worse  than  Havelock's  approximation 
over  a  large  portion  of  the  surface  of  the  spheroid.  The  approximation 
obtained  keeping  twenty. equations  in  (44)  is  depicted  in  Fig.  4,  where 
the  relative  density  is  plotted  along  meridian  planes  of  the  spheroid, 
together  with  Havelock's  approximation. 

The  solution  obtained  in  this  work  is  an  "exact"  solution  within 
the  theory  of  infinitesimal  waves.  The  determination,  for  the  same 
particular  shape,  of  the  error  associated  with  the  linearization  of 
the  free  surface  boundary  condition,  using  the  approximation  techniques 
leading  to  (25)  for  the  second-order  contribution  to  the  solution  of  the 
nonlinear  problem,  constitutes  a  natural  continuation  of  the  present 
study.  For  the  circular  cylinder,  Tuck  (1966)  has  shown  that  the  error 
associated  with  the  failure  to  satisfy  the  free-surface  boundary  condi¬ 
tion  is  larger  than  that  associated  with  the  failure  to  satisfy  the 
boundary  condition  on  the  surface  of  the  body.  The  circular  cylinder, 
however,  is  a  shape  of  little  practical  interest,  and  the  corresponding 
result  for  slender  prolate  spheroids  should  prove  a  valuable  contribution 
to  the  general  knowledge  in  this  area. 


CONCLUSIONS 


The  flow  about  a  submerged  prolate  spheroid  in  axial  horizontal 
motion  beneath  a  free  surface  has  been  treated  in  this  work  and  an 
"exact"  solution  within  the  theory  of  inf init-simal  waves  has  been 
obtained.  Comparison  of  this  solution  with  Havelock's  approximation 
reveals  that  significant  errors  are  associated  with  the  latter,  in  which 
the  boundary  condition  on  the  surface  of  the  body  is  not  satisfied 
exactly.  For  a  prolate  spheroid  with  slenderness  ratio  (the  ratio  of 
major  to  minor  axis)  slightly  larger  than  five,  a  focal  distance  twice 
the  depth  of  submergence,  and  a  Froude  number  (defined  with  respect 
to  the  distance  between  foci)  of  0.4,  the  wave  resistance  is  larger 
than  Havelock's  by  about  34%.  For  slenderness  ratios  of  3.64  and  2. 40, 
the  same  relative  depth  of  submergence,  and  the  same  Froude  number,  the 
corrections  are  as  much  as  90%  and  368%,  respectively,  of  Havelock's 
approximation  (the  spheroid  corresponding  to  the  latter  slenderness 
ratio  is  very  close  to  piercing  the  free  surface). 

The  successive  approximations  computed  in  order  to  obtain  the 
final  values  of  the  coefficients  of  the  required  expansions  in  sphe¬ 
roidal  harmonics,  and  the  corresponding  values  of  the  surface  density 
of  the  source  distribution  and  the  wave  resistance,  show  that  the 
corrections  are  small  when  only  a  few  terms  in  the  expansions  are  re¬ 
tained.  The  convergence  is  slower  for  the  surface  distribution,  which 


oscillates  about  the  solution,  and  the  approximation  obtained  by 
keeping  only  a  few  ter*  in  the  expansion  is,  for  some  portions  of  the 
spheroid,  farther  from  the  solution  than  Havelock's  approximation. 

An  infinite  set  of  equations,  essentially  equivalent  to  that  obtained 
in  this  work  for  determining  a  source  distribution  on  the  surface  of  the 
spheroid  which  satisfies  exactly  tne  boundary  condition  on  its  surface, 
was  obtained  by  Eessho  using  an  entirely  independent  derivation.  The 
coefficients  of  Bessho's  system  of  equations,  however,  appear  to  be  in¬ 
correct,  possibly  because  of  typographical  errors,  and  his  numerical 
evaluations  are  rather  inaccurate.  The  value  of  the  wave  resistance 
obtained  by  Bessho,  for  a  Froude  number  of  0.395,  a  focal  distance 
equal  to  twice  the  depth  of  submergence,  and  a  slenderness  ratio  of 
4.17,  exceeds  Havelock's  approximation  by  146%;  according  to  the  nu¬ 
merical  evaluations  reported  here,  the  correction  should  instead  be  in 
the  neighborhood  of  60  to  65%  of  Havelock's  value. 

The  solution  obtained  in  this  work  is  an  "exact"  solution  within 
the  theory  of  infinitesimal  waves.  The  determination  of  the  error 
associated  with  the  linearization  of  the  free-surface  boundary  condi¬ 
tion  constitutes  a  natural  continuation  of  the  present  study. 
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We  have  moreover  (see  Appendix  B) 


The  final  result  is  then 
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APPENDIX  B 


PROOF  OF  A  RELATION  INVOLVING  LEGENDRE  FUN  Cl  IONS 


The  relation 


(•*  x  l  *  AVI  _ Alf  *  AVI  4*1  ,  .  | 

p«^)- =H) 


required  in  Appendix  A,  is  here  derived.  S: 
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and  therefore  ' 
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constant. 


the1Srfi™ionship°btained  by  Lamb  in  hiS  Hydr°dynamics  (P-  1^2)  from 


0,(*)  =  l.-1)  ?7a)  ( 
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the  proof  of  which,  however,  is  not 
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given. 


To  evaluate  this  constant,  we  use  the  expansions 
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(Hobson,  p.  91),  and 
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(Hobson,  p.  108),  which  hold  when is  not  reel  and  between  1  and 
“1.  The  result  follows  immediately. 
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GX+8y+y2 

EXPRESSION  OF  e  a24fi2+v2  n  T 

’  ^  7  ~  0  IN  SPHEROIDAL  HARMONICS 


We 


require  the  expansion  of  the  function 
e  ,  a  -^"+72  =  o 


(1) 


ics  in  the  region  interior  to 


an 


m  a  series  of  spheroidal  harmon 

arbitrary  spheroid  centered  at  the  origin  Here  rh 

8n-  He^e  the  x  axis  is  assumed 

a  the  axis  of  the  spheroidal  coordinate 
In  the  particular  case 


system. 


this 


a  k.  P  -  Ik  cos  t,  r  =  ik  sin  t 
evasion  csn  be  obtains,  without  ^  di£ficuUy 
tionship 

* I  C — - j=  £  Ml)  £(<*©} 
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(Hobson  1931  n  aia's  j 

*  p.  414)  and  expanding 
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only  if  the  y  axis  is  the  axis  of  the  spheroidal  coordinate  system. 

The  technique  to  be  used  here  to  obtain  the  expansion  ,f  (8)  is  based 
again  on  results  of  Hobson  and  yields  a  more  general  expansion  of  (1). 

i-et  f(x,y,z)  be  a  potential  function  over  a  bounded  region  R 
containing  the  origin.  Choose  any  spheroid  "?  =  ^  contained  in  R, 
7l  b6iBg  We  obtain  the  expansion  of  f(x,y,z)  in  sphe¬ 

roidal  harmonics  by  first  expanding  the  value  of  this  function  on  the 
surface  S  of  the  spheroid  in  a  series  of  surface  spherical  harmonics, 

| ^  c  sa  «/#  j 
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v 


/  +  )  )  (C  -  ‘(J  P  6>) 
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^  /n9  r;  { 


and  then  writing 
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+  Cos^if-hB1  su*  ^v, if  \  P  '*11  /10\ 

W  IV'  , 

p-  M’"/,} 


°  /Mr  I 

(Hobson,  pp.  417  et  seqq.).  We  assume  the  series  expansion  (9)  to  be 

uniformly  and  absolutely  convergent  over  the  surface  S.  Its  coeffi¬ 
cients  are  given  by 
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It  should  be  noticed  that  the  expansion  for  the  space  internal  to  an 
ellipsoid  has  been  used-  The  expansion  for  the  external  space,  which 
involves  Q^(cosh  -^  )  instead  of  P^(cosh  *1),  cannot  be  used  since 
the  function  to  be  expanded  is  not  regular  at  infinity  and  therefore 
Hobson's  theorems  are  not  valid. 


Put 


x' 


cos  ©  , 


y'  =  sin  9  cos  ,  (12) 

z '  =  s in  9  s  in  . 

The  value  of  the  function  f(x,y,z)  on  the  surface  of  the  spheroid  is 
then  given  by 

f (ax' ,by ' ,bz ' ) 

where  use  has  been  made  of  the  relationships 


a  =  c  cosh  /r|1,  b  =  c  sinh  ^ 
and  the  repeated  integrals  in  (11)  become 
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(13) 


c3 


the  integration  being  taken  over  the  surface  of  the  unit  sphere 


,2^  .2.  ,2  . 


X  *  -ry  ~-i 


-  1 ,  and  denoting  any  of  the  2n-rl  independent 


ordinary  spherical  haraonics  of  degree  n; 


r,a  Pffl(cos  6  ),  r'B  P^Ccos  &  )  eos  a  &  ,  r,u  P^Ccos  9  >  3lri  a  Ur  f 
n  n  v  n '  "  ; 


2*1  /  n 
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with 


r'-  =  x'^-ry:“+z;' 


These  functions,  as  is  well  known,  are  homogeneous  polynomials  in 
x ’ ,  y ' ,  and  z ’ . 

To  evaluate  the  double  integral  (13)  we  use  a  theorem  of  Hobson. 
We  have  (Hobson,  p.  161) 
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X  2. 
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and,  on  the  right  sxde,  x1,  y1,  and  z’  are  put  equal  to  zero  after 


/n  f  d 


■N 
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the  operations  are  performed.  The  operator  j  is  ob¬ 
tained  replacing  x1,  y1,  and  z’  in  the  homogeneous  polynomial 


Y™(x',y',z')  by  the  differential  operators  ^T»  and 


ci 

ST’ 


64 


respectively. 

It  can  be  easily  shown  that  the  conditions  of  the  theorem  are 
satisfied  for  our  function. 


where  the  factor  i  and  the  minus  sign  within  the  braces  correspond  to 
the  n  harmonics  containing  sin  mj1,  and  the  factor  1  and  the  plus 
sign  to  the  remaining  (n+1)  harmonics,  we  obtain  finally 
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If  we  write  the  requir’d  expansion  in  the  form  (6)  the  coefficients 


are  given  by 
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^here 


.2  2  2 
A  =  a  c 


In  the  particular  case  of  interest  to  this  work,  we  have 


a  =  ik  cos  t,  P  =  k,  7  =  ik  sin  t 
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Equations  (15)  and  (16)  beccxae  then 
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where  we  have  now  put 


A  =  kc  cos  t 


Since 
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Moreover 
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Making  use  of  (21),  we  can  write  (19)  and  (20)  in  the  more  compact  form 


*'  ^  2*^"  Jc  b  &Ua&CcS<^  -r  L  U.(&-  toC&tg-'C-i-  k'C.‘^<5 iu^-j  Sua'A’J  ^ 

P ,ic9j  •f  s>i  £■  x  92  cf 

*■*  I  1 


a 


=  2TT  t 


*•»  ( 

1  w 1  1  (a«t+cwtnfei 


i  \  I  \  1 


^  j  W 

v.<- _ 


(22) 


and 


LL 


^£jCb"  /  StM  nn 


■  JScl  y 


*»)  j  '*n 

=  27 Tt  P^(^^)(|5ccU^tj_ 

l 


(sc 


I—  T  (A)  (23) 


F 

I7A 


The  coefficients  of  expansion  (6)  for  the  function 


Icy+ik(x  cos  t  +  z  sin  t) 


are  then  given  by 
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where  now 


A  =  kc 


Since 
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whfere  denotes  the  modified  Bessel  function  of  the  first  kind 

of  order  (n+%),  the  coefficients  (25)  can  also  be  expressed  in  the 
form  (7)  as  we  wanted  to  show. 
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The  factor  i  and  the  minus  sign  within  the  braces  correspond  to  the 
n  harmonics  containing  sin  mVj>;  the  factor  1  and  the  plus  sign  to 
the  remaining  (n+1)  harmonics. 

For  the  proof  of  this  formula  we  make  use  again  of  a  result  put 
forward  by  Hobson.  We  have  (Hobson,  p.  93) 
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/ 

r  'W'] 

-  /J 

i 

K 

~)l ' 

i 

2(’2/w'-t'Z  | 

\f^L^Lv  *■  -  •  •  W  *1  L  <S V- * 


(Hobson,  p.  137),  and 


1. 

>]  5iU  -  (_i)  -(”+'w) ! 

(&s6 

1  *  _  (*-«*) 

2  iv*7  !  (/W-iWij  ( 

1  2(2/wa4-2) 

M-  /WJ  - 
\ 


V  /  *  ' 

2.y(z-«imi  >•'’/)-•  "jit15 


'Y  (  2/>»t  +  2^  ^2/vmtY ) 

We  make  use  of  these  expressions  to  obtain 


S  a  3  \  aax-^hy+/bz 

n  VS’  V  5L)  e 


IS/  /vj 
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Again,  for  shortness,  put 


We  have 


f  =  Oiax+pby+ybz 


9  ?  T- 

—  e  =  «.  e 


i4-‘£ii- i)  *f-  <>vi -  w.f 

and 

/?  . ; n"+/ ?  .  J  \w)  ■?  /'*’(,  .  7 

Moreover,  since  for  all  values  of  A<  which  are  not  on  the  real  line 


segment  (-1,+1), 


p*yM  = 


/  .  \  I  ^  ,  >  (  \  — 

(/Vff^AVl  I  *  .  \  ^  (rtW/VV\  I  l  A*  -  /» —  |  I  v 

~0  /A  -  J - - - - /*  i 

Z**1  mi  '  ]/  z' 

^_/wr a | /wi 2. )  /v'~’ ^  ■*"  | 

2.V(2am  +  eJ( Z^+Vj  /*  (Z'  >  j 


(Hobson,  p.  93),  we  have 


..  nan 


EU 
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When  use  is  made  of  (6),  (7),  (8),  and  (10)  the  required  result 
follows  immediately. 


APPENDIX  E 


COMPUTER  PROGRAMS 


n  o  o  ~i  o  o 


twcww-n  mmm  asm* 


“1T3* 


■■ *  ftp  W 


ft*. 


W^VtSZZK&t^' 

f 


EVALUATION  np  Sl'iOLF  INTEGRAL  -5Y  Hfc'-V- 1 TF— GA-JSS  1 


>  ft  r 


i :  :•„ : 


DDiiBl  F  P-fCISIOM  3FSSJ,  ARC,  FR 

n  IMFMSION  X( 10) ,W( IC) ,F{  10)  ,  SINT (2000)  ,  I  VO (2.0  )-• )  , 

I  INm'-*(?OG0) 

c 

READ  {  5  » "* )  MTUAO,  |  X  « I  ),.‘(I  )  ,  1  =  1  ,':CUAD) 

3  FORMAT  (I  lO/UFV).]  5)  ) 

WRITF  (6,4)  (X(  !),',.  (I  ),  1*1,  NOUAD) 

4  FORMAT  {  IFi  1 »  1  ft  X ,  4HX  (  I  )  ,  3ft  X , 4HW{  I  )  //  {  1H  ,  r  30.  1  5 ,  I  V  ,  f 3  ) . 1  5 )  ) 
C 

PI  =  3. 141*93 
C 

READ  (5‘,7)  FRR,  LIMIT,  MUM  I  NT 
7  FORMAT  (F2Q. 10,2110) 

READ  (5,1)  FR,  OOC 

1  FORMAT  (2F10.5) 

C 

GCOU2  =  1.0/(2.0*FR**?) 

G0CU2  =  GC0U2*U0C 

WRITF  (ft,?)  FR ,  GC0U2,  GOOU2,  DOG,  F R R 

2  FORMAT  ( IH1 , 2 1HFR  =  U/ SORT { ?*G*C )  =  ,*10.5, 10X, 

1 1 1HG*C/U**2  =  , F10.5, 10X,  UHG*D/U**2  =  , F 10 . ?  ,  1 C X  , 

26H0/C.  =  ,F10.5, 10X/ 1HC  ,  6HFRR  =  ,F20.i5) 

C 

H  =  1.0/  S OR T  (  ?  .  0*G DO U 2 ) 

FACT  =  8.0*PI*  FXP(-2 .0*GDCU?)#H 
C 

DO  2°  J  = Iv MUM  I  NT 
C 

R  F AD  (5,10)  M,  M,  NP,  MP 
19  FORMAT  (4110) 

INDNUM(J)  =  1000*N+100*M+10*NP+MP 
I  NT  {  J  )  =  Nf-M+NP  +  MP 
C 

IF  ( lND{J)/2*?-IN0( J) )  15,16,15 
lft  SIMT(J)  =  0.0 
GO  TO  29 

r 

15  WRITF  ( ft , 6 1  )  N ,  M,  NP,  MP 

61  FORMAT  ( / / / 1  HO , 4HM  =  ,I2,10X,4HM  =  , I  2 , 10X , 5HMP  =  ,12, 

I I  OX , 5HMP  *  ,12) 

SUM  =0.0 


onono  o  ooooooooo 


«*  J.UMWMI 


o;i  3c  i  =  i,nq')ad 

r 

V. 

XH  =  X{ I ) *U 

PC  =  S O'-  T  (  1  .  3  +  X H* * 2 ) 

KM  =  {  PC+XH)  ***-1 
RMR  =  (RC+XH)**MP 
ARC  =  GCnil?*RC 
F  R  =  EPP 

F  (  I  )  =  C.?5*(RM+1  . P/3 ”. )  *( RMP *-1 . 0/R  VP )  £RF S3J  ( \  f  A3  ,  EF  .LIMIT) 
1*3F  SS J (  NP  ,  ARG  » ER »L I  MI T  }  /  RC 
30  SU"  =  SUM+F( I )*W( I ) 

C 

SIf  T f  J )  =  SUVSFACT 

WRITT  (6,5)  SINT(J),  SUM,  { I ,F( I ) , 1=1 ,NGUAU) 

5  FORMAT  (  1  HO ,  7 HS  I  NT  =  ,  F 20. 1 5 , 20  X ,  6HSUM  -  ,  F  7f  .  I  5/ IHC,  9X  , 
IIHI, 20X, IHF//( 1U  , I 10.F30. 153 ) 

C 

29  COf  TINUF 
C 

WRITE  (6,499)  (SINT(I),  INONUM(I),  I  =1,  N'U  •*!  NT  ) 

499  FORMAT  ( 1151, 1 7X  ,4HS  INT  ,  25X  ,4HNV.NK/// ( 1H0 , 1CX  , F2  3. 15 , 13X, 

1 1 10 ) ) 

WPITF  (7,599)  (SI  MT  (  I  ) ,  TNDNUM(I),  I  =  l,NUMINT) 

599  FORMAT  (E20.10,  20X,  110) 

CALL  EXIT 
END 


FUNCTION  5ESSJ 

BESSEL  FUNCTION  J  OF  ORDER  N+l/2 


DOUBLE  PRFCISION  FUNCTION  BESS J ( N, X , FRR , L I M I T ) 

DOUBLE  PRECISION  X,  ERR,  PI,  FACTOR,  TERM,  SUM,  X2,  DSIN, 
10C0S,  OSQRT,  DABS,  3ESJ,  FACT 
DIMENSION  BESJ(IOO) 

PI  =  3.  1415926535898 
IF  (X-5.0)  2,2,3 

3  NP1.  =  N+-1  4 

THF  RECURRENCE  RELATION 

J ( N*1 , X )  +  J ( N— 1 , X )  =  (  2*N/X  ) ( N,X  ) 

IS  USED 

IF  ( N — 1 )  13,11,16 


Iirnlr 


n 


•umiii  nmnct  trwru'wi  i*. 


o  o  o  o 


13  RESJ(I)  =  DSIN(X) 

GH  T(’  1  R 

11  RFSJ{?)  =  0S1N1  X)/X-i^nS(  X  ) 

GO  TP  14 

16  p.rsjt  n  =  isr\'(  xi 

p  £s  j  ( 2  >  =  'isiN.x)/y-"»cns(Xj 

DP  1^  1=3, NP1 
M  =  I 

19  BFSJII)  =  (2.0*AI-?.  .0  )/y*RES J<  1-1 1--3ESJ  (  1-2) 

18  RESSJ  =  P  =  SJdlPl}*OSg^T{?.«/(PI*xn 
30  TO  101 

2  FACT  =  1.0 

THf  EXPRESSION  OF  A  BESSEL  FUNCTION!  J  AS  A  S T  <  1  - S  OF 
ASCENDING  POWERS  OF  THE  ARGUMENT  IS  USED 

IF  ( N)  8,7,8 
8  DO  1  I  =  UN 
A2IP1  =  ?*I+! 

5  FACT  =  FACT*X/A2IP1 

7  FACTOR  =  DS3RT(2.0/PI)*(X**(0.5> >*FACT 
TERM  =  1.0 
SUM  =  1 .0 
AN 2  =  2*N 
X2  =  X**2 
DO  10  I  =  l, LIMIT 
AI2  =  2*1 

TERM  =  (-l.0)*TEP.M*X2/(AI.2*{  AN2+AI2  +  1  .0  )  ) 

SUM  =  SUM+TERM 

IF  ( DARS( TERMJ-DARS (EPR+SUR) )  12,12,10 
10  CONTINUE 

WRITE  (6,113)  TERM,  SUM,  X 

113  FORMAT  (  1H0,40HREQUIRE0  ACCURACY  NOT  MET  IN  LIMIT  ST!:^S/lH 
17HTERM  =  ,F20.15,IOX, 6HSUM  =  ,F 20. 1 5, 10X,4HX  =  ,F2P.}S) 

12  RESSJ  =  SUM*F ACTOR 
101  RETURN 

END 


n  o  o 


EVAIUATIPN  Cif-  SINGLE  INTEGRAL  BY  SIMPSON  S  BUI  E 

C . 

C 

LXTEEMAL  f 

enr'MPN  NS  M,  \'P,  «p,  GOPIJ?,  GCCJU2,  EXRuES,  li  mbps 
DlfTNSION  SINKIOOD,  INQUOOC),  INONUM  ( 1  303  ) 

C 

PI  =  3.141593 
C 

READ  (*,11)  ERE,  FPP.3ES,  L I  WI  Ti  ,  LIMIT?,  LI'.BSS,  Nl! "  l  NT 
11  FORMAT  (2F20. 10, 31 10/1 1G) 

READ  (5,1)  FR ,  DOC 

1  FORMAT  ( ?F 1 0. 5 ) 

C 

GCrU?  =  1 .0/( 2.0*FR**2) 

GDCU2  =  GC0U2-03C 

WRITE  (6,2)  FB%  GCPU2,  G00U2,  DOC,  ERR 

2  FORMAT  (1H1,21HFR  =  U/SOR  T(  ?*G*C )  =  ,F10.5,10X, 

11 IHG*C/U**2  =  ,F10.5, 10X, 1 1HG*D /U** 2  =  ,F13.5,10X, 

26HD/C  =  ,F10.S, 10X/1MC,6HERR  =  ,F20.15) 

C 

FACT  =  H.O*PI*  EXP  (-2  .0*0001)2 ) 

H  =  1.0/  SORT ( 2 . 0*GD0U2 ) 

C 

DO  7  1*1, NUMI NT 
C 

REAR  (5,21)  N,  M,  NP,  MP 
21  FORMAT  (4110) 

INDNUMII)  *  1000*N*100*M+10*NP+MP 
IND(l)  *  N+M+NP+MP 
C 

IF  ( I NP ( I )/2*2-!ND( I) )  25,26,25 
26  S  INT (  I)  =  0.0 

GO  TO  7 

C 

25  WRITE  (6,61)  N,  M,  NP,  MP 

61  FOPMAT  ( / / /I  HO , 4HN  =  ,I2,10X,4HM  =  ,!2,10X,5HNP  =  ,12, 

1 10X , 5HMP  =  ,12) 

A  =  0.0 

B  =  SQBT(  10.0)*H 
SUPS  =0.0 
ERUSFD  =  ERR 
C 

DO  14  K=1 ,L IM  IT  1 


—  ad 


nnonono 


50  CAtl  S«PSVIF,  A,  5,  ERUSEO,  LIMIT?,  SI  1  ,  S,  Nir*,  IF**) 

WRIT!  ( 6 , ?  0 1  )  Silt  S,  MUM,  IER,  ERUSEO 
201  FORMAT  (1H0.6HSI1  =  , F  ?0.  1  5 , 1  OX , 4HS  =  »F?Q.15»10X  =  , 

1 15, 1CX,  MilFR  =  ,  15/  H0,9HFPUSED  =  ,F?0.15) 

IF  (  I  f  R  —  4  )  5,  <5, 5 

6  FRltSED  =  FRUS Ef'*10.  0 
GO  TP  50 

C 

5  SUMS  =  SUMS  +  S 

IF  (  ARS(S)-  A«S(ERF*SUMS ) )  13,13,144 
144  A  =  B 

AKP 1  =  K+l 

14  B  =  SORT! AKPl*10i*H 
C 

WRITE  (A, 15) 

15  FORMAT  (  1HO,41MREQUIREO  ACCURACY  NOT  MET  IN  LIMITl  ST.-N>S) 

C 

13  SINT ( I)  =  SUMS* FACT 

WRITE  (6,91)  SINT(I),  ERUSED 
91  FORMAT  ( I  HO, 7HS I  NT  =  , F20. 1 5, 10X.9MFRUSF0  =  ,F?0.15) 

C 

7  CONTINUE 
C 

WRITE  (6,499)  (SINT(I),  INONUM(I),  1=1,  NU:,*I  NT  ) 

499  FORMAT  (  1H1  ,  1 7X  ,  4HS  I  NT  ,  25X  ,4HNMNM/// (  1H0 , 10  X ,  •=  ’ :) .  1  5  , }  OX, 

liion 

CALL  EXIT 
END 


FUNCTION  F  ( X ) 


FUNCTIO.1  F  (  X ) 

C 

COMMON  Nt  M,  NP,  ,MP,  GD0U2 ,  GCUU2 ,  ERRBFS,  LIMBES 
DOUBLE  PRECISION  BESSJ,  ER,  ARG,  DEXP,  GD0U2D,  X2P 
C 

X2  =  X**2 

RC  =  SORT ( I . 0+X2 ) 

RM  =  (RC+X)**M  * 

RMP  =  ( RC  +  X ) **MP 
ARG  =  GC0U2*RC 
N1  =  N 
N2  =  NP 
ER  =  ERRBES 
LIM  =  LIMBES 


nnnonno 


EVALUATION  OF  OOUDLF  INTEGRAL  3Y  SIMPSON  S  RULE 


EXTERNAL  G 

DIMENSION  DINT  (1000),  IND(IOOO),  INDNUM(IOOO) 

COMMON  N,  M,  NP,  MP,  GDDU2,  GC0U2,  ERP3ES ,  LIMoFS  ,  HDDER, 
1ERR0.R,  LIMOFR,  ERR,  LIMIT!,  LIMIT?,  KO,  FP,  RC  ,  NTMT8 
C 

PI  =  3.141553 
NTMTB  =  0 
C 

READ  (5,11)  EP»,  ERP.RFS,  LIMIT1,  LIMIT.?,  LIMBES,  NUMINT 
11  FORMAT  (2F20. 10,3110/1 10) 

READ  (5,47)  HODER,  ERRDFR,  LIMDEP 
47  FORMAT  ( F 10 . 5 ,F 20. 1 0, 1  10 ) 

RFAO  (5,1)  FR ,  DOC 

1  FORMAT  (2F10.5) 

C 

GC0tJ2  =  1 .0/(2.0*FR**2) 

GD0U2  =  GC0U2*00C 

WRITE  (6,2)  FR,  GC0U2,  G00U2,  DOC,  ERR 

2  FORMAT  ( lHl ,21  HER  =  U/SQRT ( 2*G*C )  =  ,F10.5,10X, 
.111HG*C/U**2  =  ,F10.5,10X,11HG*0/U**2  =  ,F10.5,10X, 

26HD/C  =  ,F10.5, 10X/1H0,6'JPR  =  ,F20.15) 

WRITE  (6,701)  LIMIT1,  LIMIT?,  ERRBES,  LIM5ES,  HODER, 
1ERRDER,  LIMDER,  NUMINT 

701  FORMAT  ( 1H0 ,9HL  IMIT1  =  ,I5»10X, 9HL I  MIT?  =  ,I5/IH0, 
19HERPBES  =  ,F20.15, 10X,9HLIMBES  =  ,  15 /I  HO  ,8  HHO DER  =  , 
2F10.5, 10X,9HERRDER  =  , F20. 15, 10X,9HLI MOER  =  ,I5/1H0, 
39HNUMI NT  =  ,15) 

C 

DO  7  1=1, NUMINT 
C 

RFAO  (5,21)  N,  M,  NP,  MP 
21  FORMAT  (4110) 

INDNUM(I)  =  1000*N+100*M+10*NP+MP 
IND( I )  *  N+M+NP+MP 
C 

IF  ( IND ( I ) /2*2— I N0( I ) )  26,25,26 
26  DINT (  I  )  =  0.0 
GO  TO  7 
C 

25  WRITE  (6,61)  N,  M,  NP,  MP 

61  F0PMAT(////1H0,4HN  =  ,I2,10X,4HM  =  ,I2,10X,5HNP  =  ,12, 

1 10X, 5HMP  =  ,  12//) 


r>  o  o  o  o  o  o 


r 


A  =  0.0 
B  =  5.0 
SUMS  =0.0 
cRl'SEO  =  10.0#EP.R 

DO  I '  L=l»LIMITl 

50  CAIL  SMPSN1 ( G »  At  B,  ERUSED,  LIMIT2,  SI1,  S,  MUM,  IER) 
WRITE  { 6, 20 L )  SI1,  S,  NUM,  IER,  ERUSED,  A,  5 
201  FORMAT  { // // 1  HO , 6HS 1 1  =  , F20 .1 5 , 10X,4HS  =  ,F20.15,10X, 
L6HNUM  =  , I5,L0X,6HIFR  =  ,  I  5/ 1  HO , 9HE RUSFD  =  ,F20.15,10X, 

2 4 HA  =  , F20 . 15, 1QX,4HS  =  ,F20.15////) 

IF  { I FR-4 )  5,6,5  * 

6  FRUSED  =  ERlJSED*10.0 
IF  { FRUSFD.GE.O. 1 )  CALL  EXIT 
GO  TO  50 

5  SUMS  =  SUMS  +  S 

IF  (  ABS(S)-  ABS ( ERR# SUMS )  )  13,13,144 
144  A  =  B  . 

14  b  =  B  +  20.0 
C 

B  B  =  R  —  20.0 
WRITE  (6,15)  BR 

15  FORMAT  ( 1H0,41HRE0UIRED  ACCURACY  NOT  MET  IN  LIMIT1  STEPS, 
1 10X , 4HB  =  ,F20. 15) 

C 

13  01  NT (I)  =  SUMS*GC0U2 

WRITE  (6,91)  OINT(I),  ERUSEO,  B,  SUMS 
91  FORMAT  (/////  /I  HO  ,7  HO  I  NT  =  ,F^0. 15, 1  OX,  9H  ERIJSEO  =  v 
1F2C.15,10X,4HB  =  , F 20, 1 0/1H0 , 7HSUMS  =  ,F20.15) 

C 

7  CONTINUE 
C 

WRITF  (6,499)  (OINT(I),  INONUM(T),  I=l,NUMINT) 

499  FOFMAT  ( 1H1 , 17X , 4H0 INT , 2 5X ,4HNMNM/// ( iHO , 10X , F20 . 15 , 10X , 
1110)  ) 

WRITE  (7,599)  (f)INT(I),  INDNUM(I),  I  =  1,NUMINT) 

599  FORMAT  (E20.10,  20X ,  110) 

C 

CALL  EXIT 
END 


1 


FUNCTION  G (U ) 


FUNCTION  G  ( U) 
EXTERNAL  FM 


o  r>  o  o  r> 


SEAL  KO,  K 

COMMON  NS  M  t  NP  ,  MP  ,  GOOU?,  GCOU2»  ERR3ES,  LTM8ES,  HOOFp., 
1ERPOFR,  L I MOERt  ERR,  LIMIT1,  LIMIT2,  KO ,  FP,  RC,  MTMTB 

THE  DERIVATIVE  OF  THE  FUNCTION  AT  KO  =  1.0+U**2  IS 
EVALUATED  BY  DIFFERENTIATING  AN  INTERPOLATION  POLYNOMIAL 
OF  THE  FOURTH  DEGREE 

KO  =  1.0+U**2 
RC  =  SORT ( KO ) 

EROEP  =  FRROER 

3  H  =  HODER 
PFP  =  0.0 

DO  10  J  =  l , L  IMDER 
H  =  H/2.0 

FP  =  (F(K0-2.O*H)-8.0*F{KO-H)+8.O*F(KO+H)-F(K0+2.O*H) ) 

1/(12. 0#H ) 

IF  (  A3  S(  FP-PFP  )—  ABSl FP*ERDER)  )  4,<t,l0 
10  PFP  =  FP 

WRITF  (6,60)  li,  FP,  H 

60  FORMAT  ( // // 1H0 , 41  HREQU IR ED  ACCURACY  NOT  MET  IN  L IMDER  STEPS, 
15X  ,  4HU  =  ,F20.15>-5X,5HFP  =  ,F20 . 1 5 , 5X  ,4HH  =  ,F20.15) 

if  (  abs(fp).le.i.oe-k)  go  in  4 
FRDFR  =  ERDERfrlO.O 
IF  (ERDER.GT.0.1)  CALL  EXIT 
GO  TO  3 

4  SUMS  =  0.0 
ERUSFD  =  ERR 
A  =  0.0 

B  =  0.0 
U2  =  U#*? 

00  90  I  =  l,LIMITl 
A  =  B 

B  =  B  +  2 . 5/GDOU2 
IF  ( B.GF.U2)  B  =U2 

91  CALL  SMPSN  (FM,  A,  B,  FRUSED,  LIMIT2,  SI1,  S,  NUM,  TER) 

IF  (IER-4)  77,88,77 
88  FRUSED  =  ERUSFD*10.0 

IF  ( ERUS ED.GE.O .1 )  CALL  EXIT 
GO  TO  91 

77  SUMS  =  SUMS  +  S 
IF  (B-U2)  96,92,92 

96  IF  (  A3  S  (  S  )  -  A8S(ERR*SUMS) )  97,97,90 
90  CONTINUE 

WRITE  (6,93)  U 

93  FORMAT  (IH0,32HKMAX  GREATER  THAN  2*L I  MI T 1 /G00U2  , 1  OX , 


14HU  =  t F20. i 5) 

CALL  EXIT 
C 

97  A  =  U2 

R  =  2.0  +  U2 
FRUSED  =  ERR 
[MAX  =  5 

CALL  SMPSN  [FM»  A,  B»  ERUSED,  [MAX,  Sll,  PVAL,  NIJM,  IER) 

IF  {  ABS(PVAL)-  ABS { ERR*SUMS ) )  13,92.92 
C 

92  A  =  tJ2 

R  =  2.0  +  U? 

ERUSED  =  ERR 

94  CALL  SMPSN  (EM,  A,  B,  ERUSED,  LIMIT?,  Sll,  °VAL ,  NtJM,  IER) 
*  IF  (IER-4)  55,66,55 

66  IF  (  ABS { PVAL  )—  ABS ( ERR *SUMS ) )  55,55,606 
606  ERUSED  =  ERUSFD*10 . 0 

IF  ( ERUSEO.GE.O. 1)  CALLFXI T 
GO  TP  94 

55  SUMS  =  SUMS  +  PVAL 
C 

A  =  B 

B  =  B  +  2.0/GOOU? 

ERUSED  =  ERR 
DO  14  L= l *L IM  IT  1 

50  CALL  SMPSN  (FM,  A,  R,  ERUSED,  LIMIT2,  Sll,  S,  MUM,  IER) 

IF  (IER-4)  5,6,5 

6  IF  (  ABS ( S ) —  ABS ( ERR*SUMS ) )  5,5,6666 
6666  ERUSED  =  FRUSED*10.0 

IF  (FRUSED.GE.0.1 1  CALL  EXIT 
GO  TO  50 

5  SUMS  =  SUMS  +  S 

IF  (  ABS(S)-  ABS(ERR*SUMS))  13,13,144 
1^4  A  =  R 

14  R  =  B  +  2 . 0/GDDU2 
C 

WRITE  (6,15)  U,  R 

15  FORMAT  ( 1H0,41HREQUIRED  ACCURACY  NOT  MET  IN  LIMIT1  STEPS, 

1 10X ,4HU  =  ,F20. 15, 10X ,4HB  =  ,F20.15) 

NTMTR  =  NTMTB  +  1 
IF  (NTMTR. GT. 50)  CALL  EXIT 
C 

13  RM  =  ( L.O+U/RC)**M 
RMP  =  ( 1 . O+U/RC ) **MP 

G  =  (RM+l.O/(RM*(KO**M)  )  )*(RMP  +  1.0/(RMP*('<0*’:‘M!>)  )  )/  . 

1  { RC.**  (  N+NP+2-M-MP )  )  *SUMS 
WRITE  (6,301)  SUMS,  U,  G,  B 

301  FORMAT  (1H0,7HSUMS  =  , F 20 . 1 5 , 1  OX ,4  HU  =  ,F20.15,IOX, 

1 4HG  =  ,F20. 15, 10X, 7HKMAX  =  ,F20.10) 


o  n  o  r>  r>  n  o 


RETURN 

r\c 


C  .  .  .  .........  . . . .  i 

c 

C  FUNCTION  F  K  (  K  ) 


FACTION  F*MK) 

SEAL  KOt  K 

COMMON  N»  M,  NP,  V.P,  GOOU?»  GCUU2,  FRR3ES,  LI  -^LS,  HOOER, 
IFRROEK,  L I MJE  R»  F*P,  HMIT1,  L-IMIT2,  KO ,  FP,  RC 

IF  (K-KO)  3,2,3 
*3  FM  =  F(K)/IK-KC) 

return 

2  FM  =  FP 
RETURN 
END 


FUNCTION  F»K) 


FUNCTION  F ( K ) 

REAL  KO,  K 

COMMON  N,  P,  NP,  P,P,  GDCU2 ,  GCOU2 ,  FRR3FS ,  LIMBES ,  HOOER, 
1  ERR OER,  LINGER,  £RR,  LIMIT1,  LIMIT2,  KO ,  FP,  RC 
DOUBLE  PRECISION  BESS JM,  ER,  ARG,  DEXP,  GD3U2D,  KO 

ARG  =  GCOU2*K/RC 
N 1  =  N 
,\?  =  \P 
ER  =  ERRBES 
LIP  =  L'lPBES 
GOCU2D  =  GDCU2 
KD  =  K 

F  =  (K+KO}*DEXP(-2.0-*KD*GDOU2D)*(  (GCOU2*K)**(N+NP)  )* 

1 B  E S  S  J  M (Nl, ARG, ER, LIP)*BESSJM(N2, ARG,  ER, LI M ) 

RETURN 

END 


ooono  noooo 
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c 

c . *. 

c 

C  FUNCTION  RFSSJM 

c 

C  BFSSFL  FUNCTION  J  OF  ORDER  N+l/2  DIVIDED  3Y  THE 

C  ARGUMENT  RAISED  TO  THE  N  +  l/2  POWER 

C 

C . 

c 

DOUBLE  PRECISION  FUNCTION  RFSSJMIN,X,PRR,LIMIT) 

DOUBLE  PRECISION  *,  ERR,  PI,  FACTOR ,  TERN,  SUM,  X?,  DSIN, 
IDCOS,  DSORT ,  DABS,  BESJ,  FACT 
DIMENSION  RFSJIIOO) 

PI  =  3. 14 159265 3 5 39 B 
IF  (X-5.0)  2,2,3 
C 

3  NPl  =  N+l 

THF  RECURRENCE  RELATION 

J ( N+l , X )+ J ( N— 1 , X )  =  I?*N/X)*J(N,X) 

IS  USED 

IF  ( N— 1 )  13,11,16 
13  RES  J ( 1)  =  DSINI X) 

GO  TO  18 

11  3F  S J I  2 )  =  DSINI X)/X-OCOS(X ) 

GO  TP  18 

16  3FSJI1)  =  DSIN(X) 

BFSJI2)  =  DSINI X)/X-DCOS(X) 

DO  1  "I  1  =  3, NPl 
A I  =  I 

19  RESJII)  =  I  2.0*AI-3.0)/X*9FSJI  1-1  )-RFSJ  I  1-2  ) 

18  BESS JM  =  BESJ(NPl)*0SQRT(2.0/PI}/{X**(iN  +  l)  ) 

GO  TO  101 

2  FACT  =  1.0 

THF  EXPRESSION  OF  A  OESSFL  FUNCTION  J  AS  A  SFRIES  OF 
ASCENDING  POWERS  OF  THE  ARGUMENT  IS  USE.D 

IF  IN)  8,7,8  • 

8  DO  5  I  =  1 ,  N 
A  2 1  P 1  =  2*1  +  1 
5  FACT  =  FACT/A2IPI 
7  FACTOR  =  DS0RTI2.0/PI )*FACT 
TERM  =  1.0 
SUM  =  1.0 
AN?  =  2*N 
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x?  =  X**? 

DO  10  1=1, LIMIT 
AI?  =  2*1 

TERM  =  (-1.0)*TERM*X2/(  AI2*t  AN2+AI2+1.0) ) 

SUM  =  SUM+TERM 

IF  IDABSCTERM)-DABS(FRR*SUMn  12,12,10 
10  CONTINUE 

WRITE  16,1131  TERM,  SUM,  X 

113  FORMAT  I 1H0,40HREQUIRED  ACCURACY  NOT  MET  IN  LIMIT  STEPS/1H  , 
17HTEPM  =  ,F20.15,10X,6HSUN  =  ,F20. 1 5, 10X, 4HX  =  .F20.15) 

12  BESSJM  *  SUM*FACTOR 
101  RETURN 
END 


onnnonoo.-* 


SOLUTION  OF  THE  SYSTEM  OF  EQUATIONS 


I  I -R )  A  =  R 


DIMENSION  AINTC  20*20)  «  0120,201.  INO(?0,20),  INDNUM(?0,20) 
DIMENSION  NI20.20),  M(20,20),  NP(20,20),  MP(20,?0) 
DIMENSION  fl(  20,20) ,  BB(20,20),  BC0U400),  RC20) 

C 

R EAO  (5,1  )  EPS 

1  FORMAT  (F20.10) 

C 

RE AO  (5,2)  NUMAOC,  NUHSYS,  MNUMEQ 

2  FORMAT  (3110) 

C 

READ  (5,61)  ( ( A INT ( I,J),N(I,J),M( I , J ) ,NP( I,J),MP(I,J), 
1J=1,  I),  1  =  1,  MNUMEQ) 

61  FOPMAT  ( E20.10.26X.4I 1 ) 

C 

rfRITF  (6,991)  ( (AINT( I , J ) ,N( I,J),M( I , J ) , NP(  I ,J) ,MP( I , J ) , 
IJ=1, I), 1=1, MNUMEQ) 

991  FOPMAT  I lHl/( 1H0.F20* 10 , 1 OX, 4 II ) ) 

C 

MNFQM1  =  MNUMEQ-1 
00  10  1=1, MNEQM1 

C 

IP1  =  1*1 

C 

00  10  J=IPl, MNUMEQ 
C 

A I NT ( I , J)  =  AINTIJ, I ) 

N  (  I  ,  J )  =  NP ( J  « I ) 

M(I,J)  =  NP( J , I ) 

NP ( I , J )  =  N( J , I  ) 

MP ( I , J )  =  M( J , I ) 

C 

10  CONTINUE 
C 

DO  20  1  =  1, MNUMEQ 
C 

00  20  J=l, MNUMEQ 
C 

IN0(  I ,  J  )  =  N(  I , J ) +M(  I , J ) +NP (  I»J)+MP(I  ,  J  ) 

INDNUMI I , J)  =  1000*N(I,J)+100*M(1,J)+10»NPI I , J ) +MP ( I  ,J  ) 

C 


IF  1 INOC I  ,  J)/2*2-IND( I,J) )  25,26,25 
C 

26  Oil.J)  =  AINT( I, j)**-l)**( IND* I,J)/2*NP(  IfJ>> 

GO  TO  20 
C 

25  OU,J>  =  AINT*  I  ,  J)  *(~l)  **(  ( IND(  I ,  J  )♦  1 )  /2«-N»M  I » J)  ) 

C 

20  CONTINUE 
C  . 

WRITE  (6,662) 

662  FORMAT*  1HI,  UHMATRIX  AINT////) 

WRITE  (6,62)  ( ( AINT (I,J),J=l , MNUMEQ) ,1= 1, MNUMEQ) 

62  FORMAT  ( IH0,9E 14.6/ 1H0 , 9E 14.6/IH0,98X ,2E14. 6 ) 

C 

WRITE  (6,31) 

31  FORMAT  ( lHi ,8HMATR I X  D////) 

WRITE  (6,13)  ((0*1, J),J=l, MNUMEQ), 1  =  1, MNUMEQ) 

13  FOPMAT  ( IHO, 9E14.6/1H0,9E14.6/1H0,98X,2F14.6 ) 

C 

DO  501  K=1 , NUMAOC 
C 

READ  (5,213)  A0CM1 
213  FORMAT  (F20.I0) 

AOC  =  A0CM1+1.0 
WRITE  (6,41)  AOC  ‘ 

41  FORMAT  (1H1 ,6HA/C  *  ,F10.8) 

C 

CALL  MATRI X ( 0, B B, QOOT 1 , MNUMEQ, AOC Ml ) 

C 

WRITE  (6,73?  ( (BB( I, J),J=l, MNUMEQ) ,1=1, MNUMEQ) 

73  FORMAT  I////IH0,8HMATR IX  3////(lH0, 9E 14. 6) / IHO, 9E 14.6/ 
1 1 HO  ,98X,2E14.6) 

C 

00  515  1=1, MNUMEQ 
DO  515  J=l, MNUMEQ 
515  B( I , J )  =  -1.0*BB*I,J) 

C 

DO  51  1  =  1, MNUMEQ 
51  8( 1,1)  =  l ,0-BBI 1,1) 

C 

WRITE  (6,74)  *(B(  I,  J),J=1, MNUMEQ) ,1=1, MNUMEQ) 

74  FORMAT  ( //// IHO, 12HMATR IX  ( I-B ) ////( 1H0,9E 14. 6) / 1H0 , 

19 El 4. 6/1 HO ,2E14.6) 

C 

WRITE  *6,14)  QDOTl 

14  FOPMAT  ( ////IHO ,8HQD0T1  =  ,E14.6) 

Rl 1)  =  — l.O/QDOTl 

WRITE  (6,99)  R(l) 

99  FORMAT  (IH1»7HR(1)  =  ,E14.6) 


ononononnnnoo 


C 

NUMEO  =0 

00  501  KK= 1,  NUMSYS 
C 

NUMEO  =  NUMEO  ♦  IKK+U 

C 

00  22??  1=1, NUMEQ 

2222  RCI)  =  0.0 

Rll)  =  -1.0/300T1 
C 

00  97  J=1 .NUMEQ 
00  97  1=1, NUMEO  * 

L  =  I J-l ) *NUMF  0+1 
97  BCOLIL)  =  811, JJ 
C 

CALL  GELGIR, 8C0L, NUMEO, 1 .EPS, IER) 

C 

WRITE  (6,16)  IER 

16  FORMAT  I//////1H0,6HIER  =  ,15) 

WRITE  16,17)  (RII),I=1, NUMEQ) 

17  FORMAT  I//1H0  »6HA10  =  , E14.6,5X,6HA1 1  =  , E 1 4.6/ IHO, 6HA20  =  , 

1E14.6,5X, 6HA21  =  , E 14. 6, 5X.6HA22  *  .E14.6/1H0.6HA30  =  ,F14.6, 
1  5X.6HA31  =  , El  4.6,  5X.6HA32  =  ,E14.6,  5X.6HA33  =  .E14.6/1H0, 
16HA40  =  , E  14.6,  5X.6HA41  =  ,F14.6,  5X.6HA42  =  ,E14.6,  5X, 

16HA43  =  ,E14.6,  5X.6HA44  =  ,E14.6/1H0,6HA50  =  ,E14.6,  5X, 

16HA51  =  ,E 14.6,  5X.6HA52  =  ,E14.6,  5X.6HA53  =  ,E14.6,  5X, 

16HA54  =  ,E 14. 6/1 HO, 100X.6HA55  =  ,E14.6i 

C 

WRITE  17,717)  (RID,  I,  KK,  AOC,  1*1, NUMEQ) 

717  FORMAT  (E20. 10,1 10.20X, I 10,10X,F10.6) 

C 

501  CONTINUE 
C 

CALL  EXIT 
ENO 


SUBROUTINE  MATR I  XI 0, B , 0D0T1 , NUMEO, AOCMl ) 

PURPOSE 

OBTAIN  THE  AUGMENTED  MATRIX  OF  THE  SYSTEM  OF  EQUATIONS 

( 1-8)  A  =  C 


SUBROUTINE  MATRIXCD, B.QDOTl, NUM60, AOCMl ) 


r 


DIMENSION  0(20,20),  8(20,201 

DIMENSION  POL  1 10) ,  DPOLIIO),  QOUIOI,  DQOL(IO),  P(10,10), 
10P(10,I0I,  Q( 10,10) ,  00(10,10),  DPDQ{ 10 , 1 0 ) ,  DPLOOL(IO) 

DOUBLE  PRECISION  00,  POL,  DPOL,  OOL ,  DQOL*  P,  DP,  0,  DO,  XM1, 
1DPDQ,  DPLDOL 
C 

XM1  =  A0CM1 

CALL  LEGE ( 00, POL, DPOL , QOL, ~QOL , P ,DP ,0,DQ, XM I ) 

C 

NDEG  =  5 
C 

WRITE  (6,321 

32  FORMAT  ( ////1H0, l 1HFUNCT IONS  P) 

WRITE  (6,11)  ( PQL( I),(P(I,J),J=l,I),I=l ,NDEG ) 

C 

WRITE  (6,33) 

33  FORMAT  ( // 1H0, 1 1 HFUNCT IONS  Q) 

WRITE  (6,11)  (OOL(I),(Q(I,J)«J=l,I), 1=1 , NDEG) 

C 

WRITE  (6,34) 

34  FORMAT  ( //1H0.30HDERIVAT IVES  OF  THE  FUNCTIONS  P) 

WRITE  (6,11)  (0P0L(I),(0P(I,J)»J=1,I),1=1,N0EG) 

C 

WRITE  (6,35) 

35  FORMAT  ( //1H0, 30HDER I VATI VES  OF  THE  FUNCTIONS  0) 

WRITE  (6,11)  (DQOL(I),(DQ(I,J),J=l,I),I=l,NDEG) 

C 

11  FORMAT  ( 1H0,2F25 .15/lHO, 3F  25. 15/ 1H0*  4F25. 15/( 1H0, 5F  25.15) ) 

C 

DO  87  I=l,NOEG 

DPLDOL ( I )  =  OPOL ( I )/OQOL( I ) 

DO  87  J=1 , I 

87  DPDO( I , J )  =  0P( I,J)/DO(  I,J) 

C 

WRITE  (6,3  ) 

36  FORMAT  ( // IHO , 20HPD0T  OIVIDFD  BY  ODOT ) 

WRITE  (6,11)  ( DPLDOL ( I ) , ( DPDQ( I,J),JS1, I ), 1=1, NDEG) 

C 

DO  50  J=1 , NUMEQ 

8(1, J)  =  DPLDOL (1)*D(1,J)$3. 0/4.0 
B(2,J)  =  DPDQ(1,1)*0(2,J)*3. 0/2.0 
8(3, J)  =  DPLDOL (2)*D(3,J)*5.0/4.0 
0(4, J)  =  DPDQ(2,l)*D(4,J)*5. 0/2.0 
8(5, J)  =  0P0Q(2,2>*D( 5,J)*5.0/2.0 
8(6,  J)  =  DPLDOL ( 3  )*D(6 , J  )*7 .0/4 .0 
B ( 7, J )  =  DPQQ(3,i)*D(7,J)*7. 0/2.0 
P ( 8 , J )  =  DPDQ(3,2)*D(8,J)*7. 0/2.0 
8(9, J)  ~  0PDQ(3,3>*D(9,J)*7. 0/2.0 


1 
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C 

C 

c 


B  ( 1 0 » J ) 
P(11,J) 
R(12,J) 
B( 13* J) 
BC 14, J) 
B<15, J) 
B ( 16,  J) 
Bt 17,J) 
B  { 1 8  ,  J) 
BC  19, J) 
BC20, J) 


DPLOOL (4)*Q(10,J 1*9. 0/4.0 
0P0Q(4,l)*0(U,J)*9.0/2.0 
DPDQ(4,2)*D(12,J 1*9. 0/2.0 
DPOQC  4,3 ) *0 (13, J)*9. 0/2.0 
OPOOC  4,4) *D(14,J)*9.Q/2.0 
DPLDQL (5 ! *D( 1 5, J 1*11.0/4.0 
DPDQ(5,1)*D( 16, J)* 11. 0/2.0 
OP 0Q( 5, 2 )*D( 17, J 1*11.0/2.0 
0P00(5,3)*0(18,JP.  11.0/2.0 
OPOOC 5,4) *0(1 9, J 1*11.0/2.0 
0P00(5,5 )*D( 20, J 1*11.0/2.0 


50  CONTINUE 

Q00T1  =  OOOL(l) 


RETURN 

END 


SUBROUTINE  LEGE ( 00, PO'  , DPOL, QOL , OQOL , P, DP,0 ,0Q, XM 1 ) 

PURPOSE. 

EVALUATE  THE  LEGENDRE  FUNCTIONS  AND  THEIR  FIRST  DERIVATIVES 
FOR  VALUES  OF  THE  ARGUMENT  GREATER  THAN  1  AND  DEGREES  UP  TO  5 


SUBROUTINE  LEGE ( QO , POL, OPOL , OOL ,DQOL , P,  DP, Q,  00, XM 1 ) 

C 

01  ME  NS  ION  POL  ( 10  1 ,  OPOL(IO),  QOl(IO),  DQOL(IG),  P  (  1CT,  10>  , 
10P ( 1 Q , 1 0 ) ,  0(10,10),  DQ(10,10) 

DOUBLE  PRECISION  QO,  POL,  GPOL,  QOL,  DQOL,  P,  DP,  Q,  DO, 
IX,  XMi,  X2MI,  X2M12,  X2M13,  RC,  RC3,  RC5,  RC7,  OLOG,  DSQRT 
C 

X  =  1.0+XM1 

X2M1  =  XMI* (2 .0+XM1 ) 

X2M12  =  X2MI*X2M1 
X2M13  =  X2M 12*X2M1 
RC  =  CCQRT ( X2M1 ) 

RC3  =  X2M 1*RC 
RC5  =  X2M12*RC 
RC7  =  X2M13+RC 
C 

POL (13  =  X 

POL  (2i'  =  l.  5*X2M1+1  ,0 
POL ( 3i  =  2 • 5*X*X2M1 +X 


! 


POL ( 4 )  =  4.375* X 2M1 2+5.0 *X?M 1+1.0 
POL ( 51  =  X*( 7.875*X2M12+7.0*X2Mi+l .0) 

C 

DPCL(1»  =  1.0 

DPPL (21  =  3, 0*X 

OPOL (31  =  7.5*X2Ml+6.0 

DPOL (41  =  X*( 17.5*X2M1+10.0» 

DPOL ( 5 )  =  315.0/8.0*X2M12+52.5*X2M1+15.0 
C 

P(i,l)  =  RC 
P(?,l»  =  3 ,0*X*RC 
P( 2  » 2 )  =  3.0*X2M1 
P(3tl»  =  7 • 5*RC3+6  »  0*RC 
P(3,2»  =  15.0*X*X2M1 
P(3,3)  =  15.0*RC3 
C 

P ( 4  » 1 )  =  X*( 17.5*RC3+10.0*RC) 

P(4t 2)  =  7,5* ( 7 . 0*X2M 12+6.0*X2M1 ) 

P(4,3»  =  105.0*X*RC3 
P ( 4  » 4 )  =  105.0*X2M12 
C 

P(5» 1)  =  315. 0/8.0*RC5+52.5*RC3+15.0*RC 
P(5,2»  =  52.5*X*(3.0*X2M12+2.0*X2H1) 

P(5,3)  =  52.5*19.0* RC5 +8 .0*RC3) 

P(%4»  =  945.0*X*X2M12 
P(5,5)  =  945.0*RC5 
C 

DP ( 1 ? 1 1  *  X/RC 

DP ( 2  *  1 )  =  6.0*RC+3.0/RC 

DP ( 2  1 2)  =  6. 0*X 

DP(3,1»  =  1 • 5*X* ( 15, 0*RC+4. 0/RC ) 

DP(3» 2 )  =  45,0*X2Ml+30.0 
DP(3,3)  =  45 • 0*X*RC 
C 

DP ( 4 1 1 )  =  70, 0*RC3+72 . 5*RC+1 0.0/RC 
DP(4» 2)  =  30.0*X*(7. 0*X2M 1+3.0) 

DP ( 4  »3 )  =  105.0*(4.0*RC3+3.0*RC» 

DP(4,4»  =  420.0*X*X2M1 
C 

DP ( 5  ♦  1 )  =  15.0/ 8.0* X*( 105.0 *RC3+84. 0*RC  +8 .O/RC ) 
DP ( 5  t  2 )  =  52. 5* ( 1 5. 0* X2H 12  +  1 8. 0*X2M 1  +  4.0 ) 
DP(5,3)  =  52.5*X*(45.0*RC3+24.0*RC ) 

DP( 5t 4)  =  945.0*(5.0*X2M12+4.0*X2M1) 

DP( 5? 5)  =  4725 .0*X*RC3 
C 

00  =  0.5*( DLOG( 2 .O+XM 1  )-DLOG(  XM1 ) ) 

C 

QOL(l)  =  X*QO— 1 • 0 

QOL ( 2 )  *  ( i.5*X2Ml+1.0)*Q0-1.5*X 


‘ 
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QOl.  (  3  )  =  X*( 2.5*X2M1+1. 0)*Q0-2.  5*X2M1~ l l . 0D0/6. 000  1 

m(4)  =  (4.375*X2M12+5.0*X2M1+1.0)*Q0-X*{4.375*X2M1+  j 

125.000/12.000)  1  j 

C  PI  ( 5 )  =  X*(7.875*X2m12+7.0*X?M1+1.0)*Q0-7.875*X2K12- 
19. fc25*X2Ml-137. 000/60. ODO  , 

DOOLd)  =  Q0-X/X2M1  ■.  ‘ 

Don  (2)  =  3.0*X*QO-3.0-l.0/X2f'U 

0001.(3)  =  (7.5*X2Ml+6.0  )  *00-7 . 5 *X-X/X2M l 

OQOL ( 4)  =  2.5*X*(7.0*X2M1+4.0)*00-17.5*X2M1-95.000/6.000- 
11 . 0/X2M1 

OQOL(5)  =  15. 0/8.0*(21.0*X2Ml2  +  28.0*X2M  1+8.0)  *00-10  5.  0/8.0* 
lX*t  3.0*X2Mi+2.0)-X/X?Ml 

0(1,1)  =  Q0*RC-X/RC 

Q( 2 , 1 )  =  3.0*  X*RC*00-3. 0*RC- 1 . 0/RC 

0(2,2)  =  3.0*  X231*Q0-(3.0*X**3-5.0*X)/X2M1 

0  (  3  v.  1 )  =  (  7.5*X2M1+6.0)*RC*00-7.5*X*RC-X/RC 

0(3,2)  =  15.0*X*Q0*X2Ml-15.0*X2Ml-5.0+2.0/X2Ml 

0(3,3)  =  15.0*RC3*Q0-X*( 15.0*RC-10. O/RC +8. O/RC 3) 

; 

i 

0(4,1)  =  X*(  l7.5*X2Ml  +  10. 0)*RC*Q0-17.5*RC3-95. 090/6. ODO* 

IRC- 1 .O/RC 

Q (4, 2 )  =  ( 52 . 5*X2M1 2+45 .0*X2M1 ) *00-52 .5 *X *X  2M 1-  10 .0*X+ 

12 . 0*X/ X2M  1 

0(4,3)  =  1  05. 0*X  *RC3*Q0- l05.0*RC3-35.0*RC  +  14,0/RC-8.0/RC3 
0(4,4)  =  1 05. 0*X2M1 2*Q0-X* ( 105. 0*X2  Ml -70 .0+56 .0/ X2M1-48 . 0/ 

1X2M12) 

0(5,1)  =  l5.0/8.0*(2l.0*X2M12+28.0*X2Ml+8.0 )*RC*QO-X*( 

1315. 0/R.O*RC3+26.25*RC  +  l. O/RC  ) 

0(5,2)  =  52.5*X*(3.0*X2M12+2.0*X2M1 )*Q0-157 ,5*( X2M 12+X2M 1 ) - 
1 14 .0+2 • 0/ X2MI 

Q ( 5 , 3 )  =  52.5*(9.0*X2M1+8.0)*RC3*Q0-472.5*X*RC3-105.0*X* 
1RC+28.0*X/RC-8.0*X/RC3 

0(5,4)  =  945. 0*X*X2M1 2*00-945. 0*X2M 12-3 15.0*X2M1+126.0- 
172.0/ X2MI+48.0/X2M12 

0(5,5)  =  945 .0*RC 5*00- X* ( 945. 0*RC3- 630. 0*R3  +  504. O/RC- 
1432 . 0/RC3+384 .0/RC5 1 

J 

OQ (1,1)  =  X/RC*Q0-1.0/RC+1.0/RC3 

00(2,1)  =  (6.0*RC+3.0/RC)*Q0-6.0*X/RC+X/RC3 

DO ( 2,2)  =  6.0*X*Q0-6.0~2.0/X2Ml-4.0/X2M12 

DQ (3 , 1 )  =  X*(22.5*RC+6,0/RC) *00-22 . 5*RC-13. 5/RC+ l . O/RC 3 

0Q(3,2)  =  ( 45.0*X2Ml  +  30. 0) *QO-X*( 45 • 0+4 .0/X  2M12 ) 

DQ( 3,3 )  =  45.0*X*RC*Q0-45  »0*RC- 15.0/RC+6. 0/RC3+24„ 0/RC5  • 

DO ( 4, 1 )  =  ( 70.0*RC3  +  72 • 5*RC+ 1 0. O/RC ,  -QO-X* ( 70 .0*RC+  j 

1 155.000/6 .ODO/RC-l . 0/RC3 ) 


0Q(4,2)  =  30.  O*X*(7.0*X2Ml  +  3.0)*Q0- 210.  0-*X2M  1-160.0- 2.0/ 
l  X  2  M 1  — 4 . 0  /  X  2  M 1 2 

00(4,3)  =  10‘>^0^‘(4.0<rX2Ml+3.0)*RC:<:QO-X*{420.0*^C  +  35.0/RC  + 
114.0/RC3-24.0/RC5) 

00(4,4)  =  420.0*X*X2ML*QO-140.0-420.0*X2MH-56.0/X2M1- 
132.0/X2M12-192.0/X2M13 

00(6,1)  =  15.0*X*(105.0/S.0*X2Ml+10.5+l.0/X2MI)*RC*Q0- 
11575.0/ 8.0*RC 3- 1785.0/8. 0*RC- 165. 0/4. 0/RC+1.0/RC3 

0Q(5,2i  =  52. 5* (  15. 0*X2W1 2 ♦1B.0*X2M  1+4.0)  *00-787.  ^*X*X?M1- 
1420. 0*X-4. 0*X/X2M12 

0Q (5,3)  =  157.5*Y*( 15 . 0*X2M 1+8 . 0 )*RC*Q0-236? . 5*RC 3-2047. 5* 
IRC-105.0/RC-12.0/RC3+24.0/RC5  *' 

DQ(5,4)  =  945.0*(5.0*X2M12+4.0*X2M1)*QO-4725.0*X*X?M1- 
1630.0*X+144.0*X/X2M12-192.0*X/X2M13 
*  0Q( 5,5)  =  4725  •0*X*RC3*QQ-/,725.0*RC3-l  575 .0*RC  +630.0/RC- 
1360. 0/RC3+240.0/RC5+1 92  0. 0/RC7 

RETURN 

END 


oooooooooo 


DENSITY  OF  SOURCE  DISTRIBUTION 


PLEASE  DRAW  WITH  BLACK  INK 

DIMENSION  POL  1(10)*  Pl(10,10),  P0L(10),  PU0,10) 
DIMENSION  A (2 0 • 6 ) ,  N(20),  M(20) 

DIMENSION  FACTOR! 20,181 ) ,  SIGMAU(181),  THETA! 181 ) , 
1THFDEG! 1R1 ) 

DIMENSION  DATA (400) ,  X(1000),  Y{1000) 

DOUBLE  PRECISION  XMl,  P0L1,  PI,  POL,  n 
CALL  TRAPS (—1 ,-l ) 

CALL  PLOTS ( DAT  A , 400 ,7HDENSITY) 

READ  (5,1)  A0CM1 

1  FORMAT  (F20.10J 
AOC  =  AOCM1+1.0 
A0C2M1  =  AOC Ml* ( A0CM1 +2 •  0 ) 

FAC  *  SORT! AUC2MI 1*12.566372 
WRITE  (6,31)  AOC 

31  FORMAT  ( 1H1 , 6HA/C  =  ,F10.8) 

READ  (5,2)  MNUMEQ,  NUMSYS,  NTHETA,  MUMFI 

2  FORMAT  (4110) 

WRITE  (6,32)  MNUMEQ,  NUMSYS,  NTHETA,  NUMFI 

32  FORMAT  ( 1H0, 9HMNUME  Q  =  , I  3 , 1  OX, 9HNUMSYS  =  ,I3,10X, 
19HNTHET A  =  ,13, 10X, 8HNUMF I  *  ,13) 

I  =  0 

NMAX  =  NUMSYS— 1 

DO  10  L= 1 , NMAX 
LP1  *  L+l 
DO  10  LL  =  1 » LP1 
I  *  1  +  1 
N  ( I )  =  L 
10  M<I)  =  LL-1 

WRITE  (6,33)  (N(I),M(I), 1=1, MNUMEQ) 

33  FORMAT  ( //// 1H0, 4H  N  M//(IH  ,212)) 

C 

XMl  =  A0CM1 

CALL  PNMXG1 ( POL , P, XMl ) 


DO  17  K= i , NTHET A 


NTHEMl  =  NTHETA-1 
ANTHFl  =  NTHEMl 
KMl  =  K-l 
AKMl  =  KMl 

THETA(K)  =  3.141593/ANTHMI*AKM1 
THFDEG(K)  =  180/NTHEM1*KM1 
X(K)  =  THEOEG(K) 

FACT  =  SORT  (  (  COS(THETA(K) ))*( AOCMi+l.O- 

i  CCS( THET A( K) } ) )*FAC 

XM1  =  COS (THE TACK) )- 1 . 0 
CALL  PNMXL1(P0L1,P1,XM1 ) 

DO  17  J=1,MNUMEQ 

MJ  =  K(J) 

NJ  =  N(J) 

IF  (MJ.EQ.O)  FACTOR  ( J*K  )  =  -1 .0*P0L l (NJ )/POl (NJ J /FACT 
SSIGN  =  (-!)**( MJ+i ) 

IF  (MJ.NE.O)  FACTOR(J,K)  =  SS IGN*P1 ( NJ, M J) / P ( N J , M J J /FACT 

17  CONTINUE 

NUN EC  =  0 

DO  277  K=i,NUMSYS 

NUMEO  =  NUMEQ+K 

READ  (5,7)  (A( I, K), 1=1, NUMEQ) 

7  FORMAT  (E20.10) 

WRITE  (6,77)  NUMEQ 

77  FORMAT  ( l HI ,40HCOEFF IC I  ENTS  ANM ,  NUMBER  OF  EQUATIONS  =  ,1?) 
WRITE  (6,777)  (N(I),M( I ) ,A(1,K) ,1=1, NUMEQ) 

777  FORMAT  ( ///( 1H0 , 1HA,2I2,4H  =  ,E20.10)) 

IF  (K.EQ.l)  NUMEQ  =  0 
277  CONTINJE 

Nl  =  NTHETA 
N2  =  Nl +1 
N  3  =  Nl  +  2 
X (N2 )  =  0.0 
X ( N3 )  =  9.0 
Y(N2 )  =  -0.10 
Y(N3)  =  0.02 

WRITE  (6,500)  X ( N2 ) ,  X(N3),  Y(N2),  Y(N3),  N2,  N3 
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500  FORMAT  {//////1H0»4F1C.5»2I101 

DO  27  KFI=l,NUMFI 

NUFIM1  =  NUMFI-1 

ANFIM  =  NUF  I  Ml 

KFIM1  =  KFI-1 

AKFIM1  =  KFIM1 

FIOFC  =  10O/NUFIMI*KFIMI 

FI  =  3.141593*AKFIM1/ANFIM1 

NUMEQ  =  0 


CAIL  PLOT (6. 0,-10. 5, -3) 

CAU  PLOT  (0.0,0.75,-31 

CALL  AXIS(0.0»0.0,1  1  HAGGLE  THETA  ,- 1 1 , 20 . 0.  C.  0.  x<  w  ■»  i  v,w,, 
CALL  AXI  S  (  0.0  ,0 . 0 , 7H0ENS  ITV,  7 , 1 0.0, 90 .0,  Y  (N2  J  Y(N3)  ?n  n* 
CALL  SYMBOLU.O, 8.5, 0.14, 5WI  =  ,0. 0,5) ’ ’ 
CALL  NUMBER ( 3.0, 8,5, 0.14, FIDEG, 0.0, I) 

DO  26  KK=1,NUM$YS 

NUMEO  =  NUMEQ+KK 


21 


DO  211  K=1 ,NTHETA 
SIGMAU(K)  =  0.0 
DO  21  J=l, NUMEO 
AMJ  =  M( J ) 


Y(K1|A=(SIGMAU(K)AU(K,+A(J,KK,*FACT0R(J,K,<:  C0S(  amjVfI) 


211  CONTINUE 


WRITE  (6,22)  NUMEQ, FI  DEG 

22  °F  ^UATIONS  =  ,I2,I0X,5HF1  = 

T  !  f  «SIGMAU(KJ,THETA(K),K=1,NTHETA) 
2221F20?10) J1HO?lOX,7HSIGMA/U»23X,5HTHETA//(iHO,F20.10,10X, 


, FlO . 5  1 


CALL  LINE  (X,Y,N1, 1,0,0) 

IF  (KK.EQ.l)  NUMEQ  =  0 
26  CONTINUE 


CALL  PLOT (20.0,0 .0  /— 3 ) 

27  CONTINUE 
CALL  EXIT 
END 


-5^ 
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EVALUAT ION  OF  THE  HAVE  RESIST ANCF 


DIMENSION  AINT(20,20),  P{20, 20),  I\'D(20,20),  I  NDNIJM  (  20 ,  ?0  ) 
DIMENSION  N ( 20 , 20) »  M(20,20),  NP(20,20),  MP(?0,?0) 
DIMENSION  NN( 20 ) *  MMI20),  A(20,6),  A4(20,20) 

DIMENSION  S( 20,20) ,  SUM(20) 

C 

READ  (5,201)  NUMAOC ,  NUMSYS,  f.NU*FQ 
201  FORMAT  (3110) 

C 

READ  (5,61)  ( ( AINT( I,J) ,N( I, J) ,M( I , J),NP( I, J ) ,MP( I, J ) , 
l  J=  1 , I ) , I =1 ,MNUKEQ) 

61  FORMAT  (E20.10,26X,4II J 
C 

WRITE  (6,991)  ((AINT(I,J),N(I,J),M(I,J),NP(I,J),VP(I,J), 
1J=1,I)  ,I  =  1,MNUMEC> 

091  FORMAT  (I. Hl/(  1H0*F20.10,10X,4I1  )  ) 

C 

MNEQM1  =  MNUMEQ-i 
DO  10  I® 1 , MNEQM1 
C 

IPi  =  I+l  / 

€  ' 

DO  10  J= I  PI ,  MNUMEQ 
C 

A I  NT ( I  ,J)  =  A I  NT ( J , I ) 

N  ( I » J  )  =  NP  ( J  ,  I  ) 

M< I , J)  =  MP { J , I  ) 

NP ( I , J )  =  N ( J ,  I  ) 

MP ( I  , J )  =  M( J , I ) 

C 

10  CONTINUE 
C 

DO  20  1  =  1, MNUME  Q 
C 

DO  20  J=l, MNUME Q 
C 

I ND ( I , J  )  =  N(I,J)+M( I,J)+NP(I,J)+MP(I,J) 

I NDNUM ( I , J )  =  1000*N( I , J ) +100*M (I , J )+10*NP( I,J)+MP( I,J) 

C 

IF  ( IND( I, J )/2*2— IND( I , J) )  25,26,25 
C 

0(1, J)  =  AINT(I,J)*(-1)**(IN0( I, J)/2*NP( I , J)+M(I , J) ) 

GO  TO  20 


26 


25  D( I,J)  =  0.0 

20  continue 

WRITE  (6,662) 

662  FOR” AT (1H1,11HMATRIX  AINT////) 

WRITE  (6,62)  (( AINT(I,J),J=1,MNUMEO),I=I,HNUKEO) 
62  FORMAT  ( 1HO,9E14.6/1HO,9E14.6/JHO,98X,2E14.6) 

WRITE  (6,31) 

31  FORMAT  ( 1H1 , 8HMATR I  X  D////) 

WRITE  (‘6,13)  (  (  DC  I  ,  J )  ,  J=l,  MNUME Q)  ,  1  =  1,  WNUME 0) 

13  FORMAT  (1H0,9E14.6/1H0,9E14.6/1H0,98X,2E14.6) 

I  =  0 

NMAX  =  NUMSYS-l 

DO  11  L  =  1 » NMAX 
LPi  =  L+l 
00  11  LL= 1 »LP 1 
I  =  I+l 
NN ( I )  =  L 
11  MM( I )  =  LL-1 

DO  5555  KAOC  =  1 , NUMAOC 
0 

READ  (5,5554)  A0CM1 
5554  FORMAT  (F10.5) 

AOC  =  1.0  +  A0CM1 
WRITE  ( 6,5553 )  AOC 
5553  FORMAT  (1H1.6HA/C  =  ,F10.5) 

NUMEO  =  0 


C 

DO  277  K  =  1 » NUMSYS 
C 

NUMEO  =  NUMEQ+K 

READ  (5,7)  (A( I,K),I=1,NUMEQ) 

7  FORMAT  ( E20.10) 

WRITE  (6,77)  NUMFQ 

77  FORMAT  ( 1H1 ,40HCOEFFIC I  ENTS  ANM,  NUMBER  OF  EQUATIONS  =  ,12) 
WRITE  (6,777)  (NNCI ) ;M’M  I) ,A( I , K) , I  =  1 , NUMEQ) 

777  FORMAT  (/// ( 1H0, 1HA ,2 12 , 4H  =  , £20 .10)) 

C 

IF  (K.EQ.l)  NUMEQ  =  0 
277  CONTINUE 
C 
C 


aem&KvcM 


NUMEC  =  0 

no  50  KK=  1  * NU'iSYS 

c 

NUMEC  =  NUMEQ+KK 
SUM(KK)  =  0.0 
C 

WRITE  16,102)  NUMEC 

102  FORMAT  ( 1HI,22HNUM8ER  OF  EQUATIONS  =  ,1 2/////IH0, 1  OX, 
14HTERM, 22X ,4HNMNM) 

C 

DO  41  K.=  1 ,  NUMEQ 
DO  41  1=1 ,NUMEQ  * 

C 

AA( I ,K )  =  A( I , KK  )*A ( K , KK ) 

IF  COII,K) .EQ.O.O)  GO  TO  41 
S ( I , K )  =  AA(I,K)+D(ItK) 

c 

WRITE  (6,101)  S(  I , K  )  ,  I NDNUMI I , K ) 

101  FORMAT  ( 1H0,E20.10,10X,I 10) 

C 

SUM(KK)  =  SUM  (KM  ♦  S(I,K) 

c 

41  CONTINUE 
C 

WRITE  (6,103)  SUM?  KK ) 

103  FORMAT  (////1H0,6HSUM  =  .E20.10) 

C 

IF  (KK.EQ.l)  NUMEQ  =  C 
50  CONTINUE 
C 

5555  CC  NT1NUE 


CALL  EXIT 
END 


